Sensitivity of the Horizontal Scale of Convective Self‐Aggregation to Sea Surface Temperature in Radiative Convective Equilibrium Experiments Using a Global Nonhydrostatic Model

We conduct radiative convection equilibrium experiments using a global nonhydrostatic model to investigate the dependence of convective self‐aggregation (CSA) on domain size and sea surface temperature (SST). We use a spherical domain with radii varying from 1 to 1/16 of Earth's radius, as well as SSTs between 290 and 310 K. A single aggregation occurs at small domain sizes, whereas multiple aggregations emerge at sufficiently large domain sizes. Domain‐averaged atmospheric temperature and humidity gradually increase with domain size until they reach convergence in cases of multiple aggregations. As domain size increases, surface wind speed increases, and the boundary layer becomes more humid. Because the radiative cooling in the free atmosphere also converges in cases of multiple aggregations, the surface evaporation must be limited due to the constraint of the energy balance. Thus, the convective region shifts from single to multiple aggregations to avoid increasing the surface wind speed. Our results show that the dependence of the horizontal scale of CSA on SST is not monotonic. This dependency is closely related to changes in the structure of the detrainment from the convective region at the melting level, resulting in enhanced radiative cooling at the top of the boundary layer near the convective region for cooler and warmer SSTs. This change in the circulation structure leads to increased surface wind speed with increasing domain size. This process affects the non‐monotonic dependence of the horizontal scale of CSA on SST.

in recent years. It is generally produced in any RCE experiment under certain conditions, described below, but details of aggregation are very different between models, as revealed by the RCE Model Intercomparison Project (RCEMIP; Wing et al., 2018;Wing et al., 2020).
Due to limitations in computer resources, the onset of CSA has been investigated in RCE experiments over only a narrow or small domain using cloud-resolving model (CRM) simulations. Wing and Emanuel (2014) noted that higher SSTs are more likely to lead to CSA, while many studies have shown that CSA occurs even at low SSTs (e.g., Abbot, 2014;sss;Holloway & Woolnough, 2016). When the domain size (expressed in this paper as a length that represents the square root of its area) is sufficiently large, such as in experiments using global climate models (Becker et al., 2017) or CRMs with channel-shaped domains (e.g., Posselt et al., 2012;Wing & Cronin, 2016), CSA occurs independently of SST. Several studies have argued that the mechanisms for initiation and maintenance of CSA are related to shallow circulation, which plays a role in transporting moist static energy upgradient from dry regions to moist regions. The shallow circulation is driven by various factors such as radiative cooling in the dry region (Muller & Held, 2012), cloud-radiation interactions (Coppin & Bony, 2015;Muller & Bony, 2015), or the virtual effect of water vapor (Yang, 2018b). Other positive feedbacks to CSA have also been identified, such as interaction between convection and surface flux (Wing & Emanuel, 2014), convective moisture feedback (Grabowski & Moncrieff, 2004), and entrainment (Tompkins & Semie, 2017).
In observations, CSA is often defined by appropriate geometric indices within a square region with sides measuring about 1,000 km, using satellite data (Stein et al., 2017;Tobin et al., 2012;Xu et al., 2019). These observational studies show that with an increase in the degree of CSA, there is an increase in aridity, a decrease in cloud cover, and an increase in radiative cooling. By contrast, some attempts have been made to quantify the degree of CSA at larger spatial scales, such as the inter-tropical convergence zone (Hohenegger & Jakob, 2020;Popp & Bony, 2019). For CSA over such large spatial scales, convective clusters that are more aggregated strengthen the contrast between the dry subtropic and moist deep convective regions. This relationship is similar to that found in the idealized RCE experiments with narrow domains.
Contrary to earlier studies using RCE experiments with limited domain size, recent studies have shown that in experiments using broader domains, CSA forms multiple regions. This domain size dependency is also related to the applicability of RCE experiments to the natural atmosphere. Based on satellite observation and reanalysis data, Jakob et al. (2019) argued that a domain with a spatial scale of 5,000 km × 5,000 km could be thought of as an RCE state. Therefore, RCE simulations must be conducted with a sufficiently large domain size (with a horizontal scale greater than 5,000 km) to avoid the artificial effect of domain size. For example, it has been suggested that the equilibrium state changes with domain size; domain-averaged temperature and moisture become warmer and more humid, respectively (Arnold & Putman, 2018). A horizontal scale of 5,000 km is almost the same as the length (6,000 km) of an elongated channel geometry proposed by RCEMIP (Wing et al., 2018). A "natural" scale of CSA can be defined as the horizontal scale of the convective region in a sufficiently broad domain in which multiple convective regions co-exist. This scale can also be obtained by varying the domain size in RCE simulations. In smaller domains, such as a square domain with sides measuring 100 km (a small domain in RCEMIP; Wing et al., 2018), CSA does not generally emerge. If the domain size becomes larger than 200 km, a single cluster is formed (Muller & Held, 2012). As the domain size increases to more than 5,000 km, multiple convective clusters are formed (Arnold & Putman, 2018). Patrizio and Randall (2019) referred to this domain size as the critical domain size, such that domains below this threshold allow the formation of only a single cluster. Multiple clusters are accommodated if a domain is larger than the critical domain size.
The critical domain size is also dependent on set-up and model. For example, some climate models have only one convective region (e.g., Reed et al., 2015), while others have multiple convective regions (e.g., Popke et al., 2013;Silvers et al., 2016). In RCEMIP general circulation model (GCM) simulations, some GCMs have single convective regions, and others have multiple. In RCEMIP global cloud-resolving model (GCRM) simulations with reduced Earth radius, the Nonhydrostatic Icosahedral Atmospheric Model (NICAM) has only one convective region (particularly at high SST), while System of Atmospheric Modeling (SAM; Khairoutdinov & Randall, 2003)-GCRM has more than one convective region. In RCE experiments with long channel-shaped domains, multiple convective clusters form along the length of the channel, with a quasi-uniform banded structure across the width of the channel (Wing et al., 2020).

10.1029/2021MS002636
3 of 20 Several studies have pointed out the importance of boundary-layer processes for the dependence of the horizontal scale of CSA on SST. Wing and Cronin (2016) argued that the remoistening time scale of the boundary layer determined the horizontal scales of CSA, while Yang (2018a) proposed that boundary layer height and buoyancy determine the horizontal scale of CSA. These studies were based on RCE experiments using domains with a long channel configuration. The wavenumber of the water vapor field along the long axis of the domain was defined as a measure of the horizontal scale of CSA. The horizontal scale decreased as SST became warmer across a range of SSTs (290-310 K; Wing & Cronin, 2016;Yang, 2018a), while it increased at higher SSTs (310-325 K; Yang, 2018a). These studies showed that the boundary layer height decreased and buoyancy between moist and dry regions increased with increasing SST. Other studies showed that the pressure gradient increased and the lower-level wind speeds increased with increasing domain size in simulations with a square domain (Arnold & Putman, 2018;Patrizio & Randall, 2019). The dependence of the horizontal scale of CSA on SST has been studied only in long channel domains, not in a global spherical domain.
In this study, we use a global spherical field to investigate the horizontal scale of CSA and its dependence on SST. We use a global nonhydrostatic model with an explicit cloud microphysics scheme and without a convective parameterization scheme to investigate this dependency for domain sizes in the range between R/16 and R/1, where R is Earth's radius. The motivation behind using the realistic size of Earth is to connect the results of the idealized RCE simulations to the natural atmosphere. A comparison by Satoh et al. (2016) of the RCE experiment and an aqua-planet experiment suggested that the horizontal scale of CSA may correspond to the natural scale of a convective cluster along the equator of the aqua-planet experiment, which implies that the RCE results apply to the natural atmosphere.
The paper is organized as follows. In Section 2, we describe the numerical model and experimental design. Section 3 presents the simulation results and the dependence of domain size on specific SST and evaluates the critical domain size and convective region. Section 4 discusses why the critical domain size exists, based on varying the domain size. We discuss the mechanism underlying the dependence of critical domain size on SST using the relationship between the circulation and energy budget in Section 5. Section 6 provides a summary and conclusions.

Model and Experimental Design
We used NICAM (Satoh et al., 2008Tomita & Satoh, 2004) for the RCE experiments over a sphere. We examined the dependence on spherical domain size using domains of radii 1/16, 1/8, 1/4, 1/2, and 1 of Earth; these are referred to as the R/16, R/8, R/4, R/2, and R/1 domains (equivalent to square domains with sides of length 1,411 km, 2,822 km, 5,645 km, 11,290 km, and 22,581 km), respectively. Reed and Medeiros (2016) studied RCE using various radii with fixed grid points to investigate the resolution dependency. Here, we fixed a quasi-uniform 14-km horizontal resolution by changing grid points for each radius. The number of grid points is 10 × 4 5 , 10 × 4 6 , 10 × 4 7 , 10 × 4 8 and 10 × 4 9 , from the smallest to the largest domain, respectively. This mesh size is coarser than that typically used for cloud-resolving models covering smaller domains but is widely used for NICAM global simulations . The coarser simulation reproduces stronger updrafts that mix less with the environment than those in O(km) simulations. That may encourage aggregation, but largescale organized convective structures such as Madden-Julian oscillations are generally well-reproduced without using cumulus parameterization, even at coarser horizontal resolutions, as shown by Holloway et al. (2013). We used a vertical grid with 38 levels extending to 40 km above sea level; the layer thickness gradually increased with altitude. Cloud microphysics was computed using the single-moment bulk microphysics scheme (Roh & Satoh, 2014;Tomita, 2008). No convective parameterization was used. The turbulent closure was calculated using level 2 of the Mellor-Yamada-Nakanishi-Niino model (Nakanishi & Niino, 2004). The MstrnX scheme was used for the radiative transfer (Sekiguchi & Nakajima, 2008).
Numerical simulations were conducted in the RCE configuration with uniform SST over the spherical domain without rotation. Globally constant SSTs of 290, 295, 300, 305, and 310 K were used as the surface boundary conditions. We refer to these SST experiments as SST290, SST295, SST300, SST305, and SST310, respectively. The solar insolation was fixed at 551.58 W m −2 with a fixed zenith angle of 42.05°, which is the same as those used in RCEMIP. Initial temperature and water vapor profiles were obtained after a spin-up with a single-column version of NICAM using the same SST setting. A random noise perturbation was added to the initial temperature 4 of 20 field at the lowest five layers to break the symmetry. Each simulation was run for 150 days. We assumed that all experiments reached an equilibrium state after 100 days. We used the final 50 days for RCE analysis. Figure 1 shows the column water vapor (CWV) and precipitation distributions on the last day of the SST300 experiment. CSA was realized from R/16 to R/4 as a convective cell covering connected areas of high CWV that formed almost a single region. For domain sizes larger than R/4, multiple convective regions appear. In all cases, precipitation cells have a finer structure than that of CWV regions, and multiple regions of intense precipitation exist in a convective region defined by CWV. We focus on the structure based on CWV, which has the largest horizontal scale. The critical domain size is defined as the largest domain size that allows only a single convective region. The detailed estimation is provided in Section 3.2; we can see that the critical domain size is between R/4 and R/2 in SST300 ( Figure 1). The area of the R/4 experiments corresponds almost to that of a square domain with sides of 5,000 km; thus, the critical domain size captured by the present experiment is close to that of previous studies using a square domain (Arnold & Putman, 2018;Patrizio & Randall, 2019).

Overview of the Equilibrium Solutions of the Simulations
The results of the other SSTs are shown in Supplementary Figures S1-S4. In all experiments except SST290, the critical domain size is between R/4 and R/2. In SST290, there are multiple convective regions at R/4 and a single convective region at R/8. Note that the multiple convective regions at R/4 in SST295 are temporary and become one cluster after several days. The convective region was circular for R/16 and R/8, relatively band-shaped for R/4, and a mix of circular and band-shaped for R/2 and R/1. This difference in the shape of convective clusters between domain sizes suggests that convective regions tend to be more band-shaped as domain size increases. Convective regions are not stationary but rather change their shape, and internal precipitation systems also evolve over time. The dependence of CSA on cluster shape has been discussed in previous studies, which showed the complexity of the dependence on shapes (e.g., Holloway & Woolnough, 2016). Although it may be better to consider the shape of the convective region when referring to the horizontal scale, in this study we define the horizontal scale of the convective region as the square root of the area of the convective region. Figure 2 shows the time evolution of the global average CWV in SST300. Since the initial profile is relatively moist, the domain-averaged CWV decreases at the onset of CSA (day 20) when the contrast between moist convective regions and dry subsidence regions becomes distinct. CWV then gradually increases for domain sizes of R/4 and larger. After day 100, CWV becomes almost quasi-equilibrated and increases with domain size. For R/16, CWV oscillates after day 40 with a period of about 20 days; such oscillations are known to exist in RCE results with other models (e.g., Patrizio & Randall, 2019) and may be an interesting aspect of RCE. However, we do not examine the mechanism of these oscillations in this study; we focus only on the quasi-stationary fields. We analyzed days 100-150 and confirmed that the quasi-steady state is reached after day 100 at other SSTs (not shown).
Next, we investigate the dependence of the vertical structure of global average fields on domain size. The global average temperature gradually increases in the free troposphere as domain size increases. However, this increase stops almost completely at R/2 and R/1, indicating convergence of the global average physical quantities. As the temperature increases, the tropopause rises. These temperature responses indicate that the temperature profile becomes closer to a warmer moist adiabat at a sufficiently large domain size. The relative humidity increases in the low and middle troposphere. As the tropopause rises, it becomes drier around 10-13 km ( Figure 3b). Figure 3c shows that the melting level at around 4-5 km increases with increasing domain size. The amount of condensate in both the liquid and solid hydrometeors increases with domain size. The altitude of the cloud top increases with the rising tropopause, and the anvil cloud fraction at around 12-13 km is not significantly sensitive to domain size ( Figure 3d). The sensitivity to domain size of features such as moistening and warming is the same as in previous studies (Arnold & Putman, 2018;Patrizio & Randall, 2019). The convergence of the equilibrium state is achieved larger than a critical domain size.

Critical Domain Size and Natural Horizontal Scale of Convective Region
The results shown in Figure 3 indicate a critical domain size above which physical quantities converge. A comparison of Figures 1 and 3 suggests that when multiple convective regions exist in the convergent regime, the domain size is larger than R/4. Thus, we can estimate the critical domain size from the convergence of the global average quantities. From the snapshots (Figure 1), we can infer that the critical domain size is between R/4 and R/2, but the precise value may depend on SST in particular. We estimate the critical domain size using the global average CWV and the mass-weighted column-averaged temperature (CAT) to represent the temperature and water vapor fields. We have chosen these variables as representative of the atmospheric field in an RCE system. In addition, the convergence of CWV standard deviation (CWVstd) and domain averaged precipitation was examined. Figure 4 shows the relationship between the global averages of CWV, CAT, CWVstd, and precipitation and the domain size for SSTs tested in this study. In all cases, these values almost converge for both R/1 and R/2. We can define the critical domain size using these values as the intersection of the converged value and the logarithmically fitted R/16, R/8, and R/4. The fitting was done for each of these values. This definition of critical domain size is sufficient for the present discussion but should be considered a first approximation. Ideally, one would search for a more accurate value by varying the domain size in finer increments. However, the combination of grid spacing and domain size is limited for the NICAM simulations over a spherical geometry.
The estimated critical domain sizes are summarized in Table 1. We used various quantities to estimate these values. All the estimated domain sizes show a similar dependence on SST. The critical domain size is largest in SST295 and SST300 and smallest in SST290. However, it is unnaturally large in SST300 when estimated using CWVstd and precipitation. The snapshots ( Figure S1) suggest the existence of multiple convective regions for R/4 in SST290, as can be inferred from Figures 4b, 4d, 4f, and 4h (cyan dots). For cases of warmer SST (SST305 and SST310), the critical domain size becomes smaller than that in SST295 and SST300. Although the critical domain sizes estimated using CWV, CAT, CWVstd, and precipitation are highly variable, the dependence of critical domain size on SST is similar. The existence of the critical domain size and   Multiple convective regions co-exist at a sufficiently large domain size, such as R/1 or R/2. In such situations, we can define the natural scale of the convective region as a representative scale of the convective region. We believe that the natural scale of the convective region is achieved in a simulation above the critical domain size. However, when there are multiple convective regions, they do not necessarily represent the natural scale of the convective region. In cases where we know what fraction of the region is convective, we can consider the natural scale of the convective region to be the size of the convective region at the critical domain size.
There are several methods to define convective regions. After comparing various definitions, we use CWV or column humidity, following Muller and Held (2012) and other studies, to target the largest scale of a convective region. We need to choose a threshold value of CWV above which the convective region is defined. In this study, the threshold of CWV is defined based on the correspondence between the vertical wind speed and the CWV, in order to take into account the dynamical characteristics and capture it as a simply connected area. We examined other definitions of the convective region; for example, it is also possible to define it using vertical velocity, a dynamic convective feature.
The threshold of CWV is defined for each experiment based on the relationship between CWV and vertical velocity. Figure 5a shows the mass-weighted average vertical velocities sorted by CWV in SST300. Moist categories of CWV correspond to regions of upward motion, while drier categories correspond to regions of downward motion. The boundary of the convective region is defined at the CWV point at which the mass-weighted average vertical velocity is zero. Figure 5b shows the relationship between the estimated CWV boundary value (circle) and the frequency distribution of CWV. The estimated CWV threshold corresponds to the saddle of the frequency distribution. This shape of the probability function of CWV is similar to that found in observational studies (Mapes et al., 2018;Masunaga & Mapes, 2020). As for the histogram itself, as domain size increases, the frequency of higher CWV values increases and lower CWV values decreases. This dependence of CWV on domain size is consistent with the results shown in Figure 2, in which the global mean CWV increases as domain size increases. Figure 5c shows the fraction of the whole domain that comprises the convective region, calculated by the CWV threshold for each SST and the domain size. This fraction is represented as a single convective region up to and including R/4 and as multiple convective regions for more than R/4. As domain size increases, the fraction of convective region tends to increase, except in SST295; it is almost constant at R/1 and R/2, indicating convergence. However, it should be noted that this is the fraction of the total area of the convective regions, not the scale of individual convective regions.
The natural scale of the convective region could be estimated by multiplying the critical domain size with the fraction of convective region (the average of R/1 and R/2; Table 1). This fraction is 11%-17% at R/2 and R/1 (Figure 5c). Therefore, the dependence of the natural horizontal scale of convection on SST is similar to that of the critical domain size. This dependence does not change for other convective region indicators such as meso-convective fraction (Muller & Held, 2012;Patrizio & Randall, 2019). The natural scale tends to be maximum in SST300 and SST295 and is smaller in SST290, SST305, and SST310. SST  Note. CWV is column water vapor. CAT is mass-weighted average column-averaged temperature. CWVstd is the standard deviation of column water vapor. We introduce another measure of the natural scale of the convective region to support our finding that the dependence of the scale of the convective region on SST is non-monotonic and reaches a maximum at an SST of around 300. We define single connected regions of CWV using a connected-component labeling algorithm to label grids in which CWV is above the threshold with four connectivity (i.e., two grids belong to the same label only if they share a common side). Then, we measure the horizontal scale of the convective region for R/2 and R/1 by labeling and measuring the size of single connected regions of CWV over the threshold. Figure 6a shows an example of the result of the connectivity of CWV (same time step as in Figure 1). Although the labeling separates large convective regions, convective regions that have not yet grown to the natural scale are temporarily separated from other convective regions or merged with larger ones. Large convective regions form a banded structure, while small regions are relatively circular. We defined the convective region of the 800 largest areas for R/1 by removing the effects of smaller areas that do not reach the natural scale of convection; this corresponds to 16 convective regions per time step, as we take measurements at 50 time steps. We implicitly assume that there are at least 16 natural convective regions for R/1 at one time step since the surface of R/1 is 16 times larger than that of R/4, which has a single region. We chose the top 200 areas for R/2 because the surface of R/1 is four times larger than that of R/2. Figure 6b shows the frequency distributions of the number of convective regions as a function of their scale, defined by the square root of the area of convective regions for R/1. There are many clusters below the scale of 500 km, as indicated by dotted lines. Solid lines in Figure 6b show the histograms of the 800 largest areas. The maximum scale of the convective region is about 6,000 km in SST305; this may be an exceptional case due to the temporary merging of convective clusters. Figure 6c shows the dependence of the mean and the median scales of the 800 (200) largest convective regions for R/1 (R/2) on SST. The scale increases from an SST of 290 K to an SST of 300 K, at which it reaches a maximum, and decreases toward an SST of 305 K. The largest mean and median scales are in SST310 in R/2; this is different from the results shown in Table 1, indicating that ambiguity remains for definitions of the scale of convective regions. Both the mean and the median scales differ between R/1 and R/2 by about 500 km, although there is a convergence of a physical quantity for domain sizes larger than the critical domain size (Figure 4). When the domain size is larger than the critical domain size, the convective region's horizontal scale is within a range that is smaller than the maximum possible scale. In addition, it is possible to increase the number of convective regions. As a result, the size of each convective region at R/1 is smaller than at R/2. However, the dependence on SST is the same except for SST310 at R/2. We examined the dependence of the critical domain size on SST using multiple methods, such as by estimating from the convergence of physical quantities or connected convective region. This dependence on SST is the same as that of the natural scale of convective region. The natural scale of convective region is between 2000 and 6,000 km for all indices. This result suggests a limit to the extent to which the convective region can expand spontaneously. The natural scale of the convective region is relatively large in SST295 and SST300 and small in SST290, SST305, and SST310, indicating that it may not increase or decrease monotonically with SST.

Why Does Critical Domain Size Exist?
We have seen that multiple convective regions co-exist under a sufficiently large domain size such as R/1 or R/2 for any SST experiment. In this section, we investigate why the critical domain size exists. In other words, we examine why a single convective region has an upper size limit and cannot expand to a large single convective region at domain size R/1 or R/2.

Spatial Variance of FMSE Budget
Following Patrizio and Randall (2019) regarding the relationship between CSA and domain size, we use the budget equation for the spatial variance of the frozen moist static energy (FMSE). The FMSE is defined by where c p is isobaric specific heat, T is temperature, g is acceleration due to gravity, z is height, q v is water vapor mixing ratio, q i is ice condensate mixing ratio, L v is latent heat of vapourization, and L f is latent heat of fusion. The column net radiative flux (RAD), the sum of latent and sensible heat release from the surface (SEF), and horizontal convergence of the flux of h are the only terms that can alter column integral of h: The angular brackets denote a density-weighted vertical integral from the surface to the top of the model. Following Wing and Emanuel (2014), we can derive Here, a prime denotes a departure from a global average value. Globally averaging and normalizing each term by ⟨ℎ⟩ ′2 yields the growth rate of 〈h〉 variance.
We focus on the quasi-equilibrium state in this study and compare the global average quantities from day 100 to day 150. Figure 7 shows the global average FMSE budget terms. Our study is not much different from previous RCE studies from the perspective of the FMSE budget. The term RAD is always positive, indicating that it contributes to the maintenance of CSA, while the convergence term is negative, contributing to damping. SEF is primarily positive but decreases as domain size increases and is negative at some SSTs or domain sizes. The difference in RAD between simulations for the same SST is slight. The atmospheric longwave cooling in the convective region decreases due to clouds, and increasing moisture relative to the dry region, due to increasing shortwave heating, tends to maintain CSA. The relative contribution of the surface flux term in the total diabatic process becomes weaker than that of the radiation term as domain size increases.
The convergence term is negative, and its magnitude approaches zero as domain size increases. The convergence term indicates that the FMSE is exported from the convective region, but it decreases as domain size increases. Patrizio and Randall (2019) hypothesized that the critical domain size might occur where the convective region switches from exporting FMSE to importing FMSE. They observed a similar decrease in FMSE export as domain size increases. However, it does not change sign, that is, we do not have a case of FMSE being imported to the convective region above the critical domain size. The decrease in FMSE export from the convective region occurs regardless of SST. The decrease in FMSE export is related to an increase in detrainment in the middle layer, as shown below. The decrease in FMSE export makes the convective region wetter when domain size increases.

The Domain Size Dependence of Boundary Layer Process and Energy Balance
We analyzed the lateral structure of the boundary layer sorted by CWV, following the analysis method used in previous studies (e.g., Bretherton et al., 2005). This analysis averaged physical quantities over a 4 × 4 grid (approximately 56 × 56 km). Here, we show only the results of SST300; the dependence of the equilibrium states on domain size is similar for the other SSTs tested (not shown). Figure 8 shows the lateral distributions of the quantities near the surface as a function of CWV. As domain size increases until R/2, the 10 m wind speeds (Figure 8a), the temperature and the specific humidity at the lowest model level (Figures 8b and 8c) and the latent heat flux (LHF; Figure 8d) increase, but they are almost unchanged at R/2 and R/1. These findings are consistent with those shown in Figure 3. The increase in 10 m wind speed until domain size reaches R/2 is due to the increase in the subsidence region when domain size increases. Warmer temperature in the convective region than in the drier region increases the horizontal pressure gradient, which needs to sustain enhanced wind speed. These tendencies are consistent with previous studies conducted on square domains (Arnold & Putman, 2018;Patrizio & Randall, 2019;Yang, 2019). Despite the moistening in the boundary layer, the LHF increases due to the increase in wind speed. The upper limit of LHF is proportional to the domain-averaged precipitation, which is constrained by the net radiative cooling. We conclude that such a tropospheric energy balance constrains the surface wind speed and determines the critical domain size. We show here how the critical domain size is determined under the constraint of the energy budget. Assuming that domain-averaged column-integrated radiative cooling is Q and surface evaporation is E, we can approximate the following balance in the RCE by neglecting the surface sensible heat flux: where A is the domain size. Sensible heat flux accounts for a relatively large fraction of the surface energy flux when SST is low (Figures 9e and 9f; e.g., SST290), but still less than latent heat flux. In addition, the change in sensible heat flux with domain size is smaller than the change in the latent heat flux. Therefore, the sensible heat flux can be neglected. The evaporation is calculated using the bulk wind speed U, where C is the bulk coefficient, and Δq is the difference in water vapor between the surface and lowest model level. We have seen that U increases with domain size (Figure 9a). Let us assume that this dependency is approximated as a linear relationship: where B is the proportionality coefficient. This relationship is also used by Arnold and Putman (2018). At the critical domain size A c , we denote Q(A C ) = Q c and Δq(A C ) = Δq c . Then critical domain size is given by:  dependence of Q and Δq on domain size, respectively; Q increases and Δq decreases with increasing domain size for all SSTs. Precipitation and latent heat fluxes also increase to balance this change (Figures 9c and 9d). Both Q and Δq increase with SST at any domain size.
The reason why the critical domain size exists can be summarized as follows. In a single convective region for domain sizes less than R/2, surface wind speed increases as domain size increases because the distance between the convective and dry regions becomes large. The change in surface wind speed matches the increases in temperature and humidity. The vertical temperature profile (Figure 3a) and the height of the tropopause increase, but they are constrained by the moist adiabat. Thus, the warmest temperature profile is constrained by moist adiabat and surface temperature. That places an upper limit on radiative cooling if relative humidity does not change significantly. The latent heat flux will also have an upper limit because of the balance with the radiative cooling (Figures 9b and 9f). The latent heat flux increases with domain size despite the decrease in Δq (Figure 9c). Therefore, latent heat flux is controlled mainly by wind speed. The upper limit to the latent heat flux and the strong relationship between latent heat flux and surface wind speed indicate that the surface wind speed must have an  Table 1, and the vertical dashed lines indicate the critical domain size obtained in Table 2. upper limit. This, in turn, indicates that domain sizes allowing only a single convective region must also have an upper limit. When domain size is greater than the critical domain size, another convective cluster needs to form, because the circulation caused by a single cluster cannot cover the whole domain. It could still be possible that a larger cluster exists at a large domain size, organized in a long band (snake-like structure), as then the distance between the convective and dry regions would still be relatively small. However, suppose the individual convection in the convective region (CWV field) behaves randomly. We speculate that the band-shaped convective region may fluctuate and merge with other convective regions at a point of large curvature, then become close to circular. Therefore, we believe there may be a temporary banded structure, but it is not stable. It may produce the critical horizontal scale of the convective region when time is averaged.

The Dependence of Critical Domain Size on SST
We found that the critical domain size does not monotonically change with SST. This section discusses how the dependence of critical domain size on SST is determined. As investigated in the previous section, the dependence of surface wind speed on domain size is essential for determining the critical domain size. We show that such a dependency is closely related to the effective structures of overturning circulation and its dependence on convective activity and SST. We constructed the mass stream function following Bretherton et al. (2005) to investigate the circulation structure. The effective mass stream function ψ is calculated as an integral over vertical velocity with respect to the CWV percentile: where ρ, w, and z are density, vertical wind speed, and height, respectively; the subscript i denotes the ith CWV bin. This method enabled us to characterize the convective detrainment associated with convective mass circulation, although we should note that the topology of the circulation was destroyed in our analysis, similar to that described by Bretherton et al. (2005).
The structure of the circulation can be determined approximately from the distribution of radiative cooling. As SST increases, the height of the melting level increases, as does the height of the detrainment of water vapor from the convective region. Mid-tropospheric detrainment from the convective region is caused by increased mid-level static stability, which in turn is due to the melting process described in tropical observational studies (e.g., Johnson et al., 1999;Yasunaga et al., 2006) and RCE studies (e.g., Patrizio & Randall, 2019;Posselt et al., 2012). The detrainment enhances the vertical moisture gradient at the melting level. Because radiative cooling is related to a vertical moisture gradient ( Figure S5), it forms a peak around the melting level and the top of the boundary layer. In addition, there is a minimum of radiative cooling below the melting level. There is a strong correlation between the height of the melting level and the height of the peak in the radiative cooling at the middle troposphere. Thus, the midlevel peak of radiative cooling shifts upwards as the melting level increases with increasing SST (Figure 10). In SST290, the melting level is close to the top of the boundary layer (about 2,000 m). In SST295 and SST300, the melting level is about 3,000-5,000 m and is detached from the top of the boundary layer. Even though the boundary layer and the melting level are separated, they have a shallow circulation in the moisture space from the melting level to the boundary layer. The height of maximum cooling at the top of the boundary layer should coincide with the height of the minimum radiative cooling just below the melting level. Hence, the cooling at the top of the boundary layer is weak near the convective region. In SST305 and SST310, the melting levels are about 7,000 m and 9,000 m, respectively. The interval between the melting levels and the top of the boundary layer increases since the boundary layer height decreases rather than increases. The radiative cooling near the melting level is no longer associated with the low-level circulation due to the rising melting level. The rising melting level leads to stronger radiative cooling at lower levels. We assume that the change in the structure of the radiative cooling in the CWV height cross-section is related to the change in the wind speed in the boundary layer. This assumption is inspired by a mass balance between the vertical wind speed and the surface flow from the dry to the convective region (Arnold & Putman, 2018;Wing & Cronin, 2016;Yang, 2018a) and radiative-driven circulation (e.g., Bony et al., 2016). The vertical advection of potential temperature is almost balanced with radiative cooling when the temperature is horizontally homogeneous under RCE simulations. The vertical wind speed in the dry region is controlled by radiative cooling and static stability. Thus, there is a relationship between radiative cooling and surface wind speed.
The radiative cooling can explain the vertical wind speed at the top of the boundary layer. The radiatively driven vertical wind speed w r can be estimated as follows (Bony et al., 2016;): where Q R and θ are radiative heating and potential temperature, respectively. Figure 11a compares w r with the actual vertical wind speed at the top of the boundary layer. We chose 1,853 m height because this height is near the top of the boundary layer or outside the boundary layer in all SST experiments. The vertical wind speed almost agrees with w r in the dry region when the CWV percentile is less than 70.
To estimate horizontal wind speed in moisture space in the boundary layer u r , we follow Wing and Cronin (2016) and assume that u r depends on the distance from the central point of the subsidence region. We regard the subsidence region as taking a circular form in which the convective region is located at the edge. From the mass continuity, u r can be calculated by: where r is the radius of the circle, and the integration is done with the moisture space (CWV height cross-section). We assumed u r (0) = 0, where r = 0 is the driest region. Figure 11b shows that u r almost agrees with the 10 m wind speed in the dry region with a CWV percentile less than 60. The 10 m wind speed is enhanced more than u r near the convective region at CWV percentiles around 80 in SST290, SST305, and SST310. Since the water vapor gradient in the lower layer does not depend significantly on SST (Figure 11c), the increased 10 m wind speed near the convective region mainly enhances the LHF anomaly (Figure 11d). We speculate that the enhanced surface wind speed near the convective region (SST290, SST305, and SST310) invigorates convective activity. As convection is strong, moisture detrainment from the convective region is enhanced near the melting level. The convective region exports less FMSE when detrainment occurs at a higher level because the middle-level humidity is lower than the boundary layer humidity. This effect increases the horizontal convergence term of Equation 3 with increasing domain size. This effect may be more pronounced in SST290, SST305, and SST310 because the wind speed is enhanced near the convective region as domain size increases in these experiments. That process may lead to a smaller critical domain size by enhancing convective activity, which increases the horizontal convergence of FMSE.
The change in the structure of surface winds affects parameter B of Equation 7, which is used for the energy balance in the previous section. The coefficient B is derived from the dependence of bulk wind speed U on area. We use Figure 9a to  Here, we explain how the surface wind speed pattern in moisture space, which is relatively enhanced in the convective region, is related to larger B. The moisture supply to the convective region increases and convective activity is enhanced when domain size increases. This process works more efficiently when the wind speed increases near the convective region. Enhancement of convective activity leads to an increase in temperature in the troposphere, contributing to an increase in the radiative cooling Q. This process implies a more significant increase in bulk wind speed to balance the increased radiative cooling. Therefore, B-that is, the rate of increase in bulk wind speed relative to increasing domain size-is larger for the circulation structure with stronger wind speeds in the convective region, such as in SST290, SST305,  Note. We calculate Q c and Δq c as the mean of R/1 and R/2. B is estimated from the fitting of Figure 10a. A c is calculated by using C = 1.0 × 10 −3 [m s −1 ] and L v = 2.5 × 10 6 [J kg −1 ]. and SST310. In general, Q/Δq decreases with increasing SST because Δq increases more, according to Clausius-Clapeyron's law, than radiative cooling related to mass flux (Satoh, 2013;Satoh & Hayashi, 1992;Vecchi & Soden, 2007). As Q/Δq decreases with increasing SST, this contributes to the decrease in critical domain size if the change in B is small.

Summary and Discussion
In this study, we investigate the formation and dependence of the critical domain size of CSA on SST in RCE experiments over a spherical geometry with various domain sizes up to Earth's radius R. We conducted RCE experiments without planetary rotation using a global nonhydrostatic model, NICAM, with various domain sizes and global uniform SST. We used NICAM with a horizontal grid size of about 14 km without cumulus parameterization to calculate deep convection explicitly. We focused on the quasi-equilibrium states and analyzed the dependence of the equilibrium states on domain size and SST.
The simulation results show a single convective region when domain size is smaller than or equal to R/4, and multiple convective regions when domain size is greater than R/4 ( Figure 1 and Figure S1). As domain size increases up to R/4, the domain-averaged humidity and the domain-averaged temperature increase. We found that the domain-averaged physical quantities converge to similar values for both R/1 and R/2. In particular, the temperature profile is constrained by a moist adiabat at the critical domain size. This convergence was confirmed for all SSTs. For domain sizes of R/1 and R/2, the global average statistics are almost the same. As domain size increases, the surface wind speed increases as long as a single convective region forms CSA, and the surface latent heat flux increases. The radiative cooling in the free troposphere is controlled by the moisture profile, which is limited by the temperature if the relative humidity does not change significantly in these RCE simulations. Because the vertical temperature profile has the warmest profile constrained by moist adiabat and surface temperature (Figure 3a), a critical domain size must exist to limit the increase in the surface wind speed related to latent heat flux (Figure 9b).
Our study found an energy constraint in addition to the momentum balance described in previous studies (Arnold & Putman, 2018;Yang, 2018a). Previous studies argued that boundary layer flow requires a horizontal pressure gradient due to temperature or moisture anomaly in the boundary layer to balance momentum loss through the surface. Arnold and Putman (2018) pointed out that the larger horizontal scale produces a stronger pressure gradient, enhancing the boundary layer wind speed. We incorporate this fact in a simplified form in our study (Equation 6). In addition, our study shows that there is an upper limit to the surface wind speed because the surface latent heat flux must be approximately balanced by the radiative cooling that occurs in the free troposphere if the surface sensible heat flux is negligible. With regard to the dependence of the horizontal scale on SST, Yang (2018a) suggested that the pressure gradient caused by the virtual effect of water vapor, which depends on SST, is important to the horizontal scale of CSA. However, the effect of other processes-such as momentum dissipation time scale or diabatic processes-on the dependence of the horizontal scale on SST is unclear. Our study suggests that a change in circulation structure with increasing SST caused by the changes in diabatic heating affects the horizontal scale of CSA.
For the dependence on SST, we found that the critical domain size does not monotonically change with SST (Table 1). It is largest at an SST of around 300 K, and smaller at cooler (290 K) and warmer SSTs (305 and 310 K). We defined various indices for convective aggregation to quantify horizontal scales of the convective region and determine the natural scale of the convective region, which is the horizontal scale of the convective region at the critical domain size. The natural horizontal scale of convection is almost proportional to the critical domain size, although the fractional area of convection changes with SST and domain size.
To further explain the existence of critical domain size, we analyze the FMSE budget (Figure 7). The increase in the convergence term of the FMSE budget as domain size increases is associated with decreasing FMSE export from the convective region, which is consistent with Patrizio and Randall (2019). Thus, there is no increase in FMSE export from the convective region with increasing domain size, as Arnold and Putman (2018) suggested. Patrizio and Randall (2019) hypothesized that the critical domain size occurs where the convective region switches from exporting FMSE to importing FMSE. However, in our simulations, FMSE export from the convective region is always positive. We argue that the dependency of FMSE export is not essential in determining the critical domain size. We hypothesize that the critical domain size occurs when the convective region becomes close to the warmest moist adiabat, where the SST limits the temperature of the lower layer, and the convective region cannot be further moistened.
We further analyzed the structure of the overturning circulation of the moisture-height space to investigate the dependence of the critical domain size on SST ( Figure 10). Water vapor export occurs near the height of the melting level. This level increases monotonically with increasing SST. However, the associated change in the pattern of radiative cooling at the top of the boundary layer in the moisture-height space is not monotonic. This radiative cooling influences the wind speed in the boundary layer and the structure of the overturning circulation in moisture space ( Figure 11). As SST becomes cooler or warmer than 300 K, a lower circulation is strengthened, and the surface wind speed becomes stronger near the convective region. These changes are fundamentally caused by the melting level being high enough when SST is higher than 300 K. If the melting level is high enough, the vertical moisture gradient is enhanced in the lower levels, which results in strong radiative cooling at lower levels and higher surface wind speed near the convective region. This surface wind pattern near the convective region enhances the dependence of surface wind speed on domain size and leads to a smaller critical domain size. We found that the pattern of circulation changes to a pattern characterized by stronger surface winds near the convective region when SST increases from 300 to 305 K. We speculate that such changes in the pattern of circulation depend on the model, although we think that the non-monotonic dependence on SST exists in all numerical models. This SST dependency would be strongly affected by the model's cloud microphysics and radiation schemes. Wing and Cronin (2016) confirmed by changing radiation schemes that the horizontal scale of CSA is modified. However, we emphasize that the circulation structure in moisture-height space is important. In this study, we show that the structure of the radiative cooling, which is affected by the detrainment height associated with increasing SST, may control the circulation structure and affect the horizontal scale of CSA.
This study clarifies the dependence of the horizontal scale of CSA on SST. The results show the non-monotonic dependence on SST, which is different from the results of previous studies (e.g., Wing & Cronin, 2016;Yang, 2018a). Of course, differences in the numerical models, such as resolution and physical parameterization, may affect the dependence of the critical domain size on SST. We speculate that differences in the domain shape of the models also have an effect. Previous studies used a domain with a channel shape, which may have trapped the dynamic wave structure more easily than a sphere shape. In addition, the changes in the structure of the effective circulation in moisture space that we point out in the present study may have been smaller in previous studies. We found that this dependency is closely related to the structure of the overturning circulation in moisture space through the increased rate of surface wind speed with increasing domain size. As domain size increases, the dry subsidence region becomes larger, so the horizontal wind from the boundary layer to the convective region increases due to the law of mass conservation. Previous studies have described this feature (e.g., Patrizio & Randall, 2019). An important factor in controlling the increased surface wind speed rate is the change in the detrainment height from the melting level and resulting radiative cooling profiles. The details of responses may depend on cloud microphysics and convective parameterization schemes. Therefore, further studies on the modulation of the horizontal scale, such as investigating the convective parameterization dependence and the response in higher resolution models, are needed.
A better theoretical understanding of the horizontal scale of self-aggregation is necessary. Recently, Yang (2021) proposed a simple model for the spatial scale of self-aggregation using a shallow water system. His model ultimately expresses spatial scales by gravity wave velocity, damping time scale, and convective density. Convective density is simply modeled by the a small-scale mass sink which is balanced by the mass source, radiative cooling. The present study shows that the effect of the radiative cooling structure on the horizontal scale of convective region. We speculate that the vertical structure of radiative cooling may control the horizontal scale by affecting the strength of convection. Further study is needed to link the simple model (e.g., Yang, 2021) with CRM simulations. In addition, the horizontal scale of CSA is so large that it includes the gradient of SST and shortwave radiation in the real world. Future work should investigate horizontal scales in such cases using experiments with interactive ocean or prescribed SST gradients.

Data Availability Statement
The data package and source code for the model used in this study are available at http://nicam.jp/. The analysis data used in the manuscript is available at https://doi.org/10.5281/zenodo.4777315, and full output is available by request. We thank three anonymous reviewers for comments that helped us substantially improve this manuscript. The research was supported by Program for Promoting Technological Development of Transportation (Ministry of Land, Infrastructure, Transport and Tourism of Japan) and by MEXT (JPMXP1020351142) under the Program for Promoting Researches on the Supercomputer Fugaku (Large Ensemble Atmospheric and Environmental Prediction for Disaster Prevention and Mitigation). This research was conducted using the Fujitsu PRIMERGY CX600M1/ CX1640M1 (Oakforest-PACS) in the Information Technology Center, The University of Tokyo.