Numerical Diffusion and Turbulent Mixing in Convective Self‐Aggregation

Convective Self‐Aggregation (CSA) is a common feature of idealized numerical simulations of the tropical atmosphere in Radiative‐Convective Equilibrium (RCE). However, at coarse grid resolution where deep convection is not fully resolved, the occurrence of this phenomenon is extremely sensitive to subgrid‐scale processes. This study examines the role of mixing and entrainment, provided by the turbulence model and the implicit numerical diffusion. The study compares the results of two models, WRF and SAM, by varying turbulence models, initial conditions, and horizontal spatial resolution. At a coarse grid resolution of 3 km, the removal of turbulent mixing prevents CSA in models with low numerical diffusivity but is preserved in models with high numerical diffusivity. When the horizontal grid resolution is refined to 1 km, CSA can only be achieved by increasing explicit turbulent mixing, even with a small amount of shallow clouds. Therefore, the sensitivity of CSA to horizontal grid resolution is not primarily caused by the decrease in shallow clouds. The analysis of the total water path spectrum suggests that the amplitude of initial humidity perturbations introduced by convection in the free troposphere is the key factor. This amplitude is regulated by turbulent mixing and diffusion at small scales. Prior to the onset of CSA, increased mixing makes updrafts more sensitive to the dryness of the free troposphere, which strengthens the moisture‐convection feedback. This leads to an increased distance between convective cores and a stronger humidity perturbation in the free troposphere, which can destabilize the RCE state.


Introduction
In the absence of lateral energy transport, the atmosphere would be in a statistical Radiative Convective Equilibrium (RCE), where radiative cooling is balanced by convective heating.Despite large-scale dynamical forcings are ubiquitous in the real atmosphere, RCE can be observed in the tropics on a daily scale, when considering areas larger than 5,000 km 2 (Jakob et al., 2019).
Without heterogeneities in boundary conditions and forcing, idealized numerical simulations of RCE provide a simple framework to study the internal interactions between moist convection, radiation and circulation.In these simulations, the atmosphere is destabilized by a combination of surface heating and tropospheric radiative cooling.The Weak Temperature Gradient (WTG) hypothesis applies to the free troposphere.Heating anomalies directly generate large-scale circulation that moistens (dries) the atmosphere through large-scale water vapor convergence (divergence).Deep convection develops uniformly throughout the domain, naturally reflecting the homogeneity of the external forcing.However, after a few days, numerous numerical studies using various Cloud-Resolving Models (CRMs) (C. S. Bretherton et al., 2005;Held et al., 1993;C. J. Muller & Held, 2012;Tompkins, 2001;Wing & Emanuel, 2014;Yanase et al., 2020) and General Circulation Models (GCMs) (Wing et al., 2020) have shown a transition from the initial uniform state of convection to a state in which convection clusters into a moist patch surrounded by a dry patch with suppressed convective activity.This process is known as Convective Self-Aggregation (CSA), since it is caused by internal feedback between moisture, convection and radiation.When CSA occurs, the atmosphere dries out, large-scale circulation develops and Outgoing Longwave Radiation (OLR) increases.These changes have significant implications for climate.Additionally, in the presence of rotation, CSA organizes into eddies that can serve as an idealized framework to study Tropical Cyclone Genesis (Carstens & Wing, 2022;Nolan et al., 2007;Ramírez Reyes & Yang, 2021;Rappin et al., 2010;Wing et al., 2016).
Idealized numerical simulations have revealed several physical processes that are relevant for the onset and maintenance of CSA (Wing et al., 2017).Radiative and surface flux feedback positively contribute to clustering in the early stages of aggregation (C. S. Bretherton et al., 2005;Tompkins & Craig, 1998).C. J. Muller and Held (2012) emphasized the importance of low-level longwave cooling in dry regions both from clear-sky and shallow clouds.This differential radiative cooling between moist and dry regions drives a low-level upgradient Moist Static Energy (MSE) circulation which expands dry regions and maintains self-aggregation, as also confirmed by Wing and Emanuel (2014).Therefore, C. J. Muller and Held (2012) suggested that reduced lowlevel cloudiness at finer horizontal resolutions may explain the absence of CSA in simulations with horizontal grid sizes below 1 km.These numerical results were used by Emanuel et al. (2014) to develop a radiativeconvective instability theory.This theory demonstrates how, in a WTG framework, an anomalous radiative cooling in the free troposphere (which corresponds to a dry perturbation) generates subsidence motions that dry the atmosphere and further increase radiative cooling there, creating a feedback loop.
In addition to radiative and surface flux feedback, other studies have demonstrated the importance of boundary layer diabatic feedback.Jeevanjee and Romps (2013) found that inhibiting cold pools (by switching off rain evaporation) allowed self-aggregation to occur at all domain sizes.C. Muller and Bony (2015) found similar results even when radiative cooling rates were homogenized across the free troposphere.This "moisture-memory aggregation" can be replicated by the simple coarsening model by Craig and Mack (2013) and has recently been found to occur spontaneously in models with a nearly saturated sub-cloud layer (Cerlini et al., 2023).Additionally, D. Yang (2018) showed that CSA can occur also with minimal ingredients by homogenizing radiative cooling throughout the entire column, switching off rain evaporation over the boundary layer and without interactive surface fluxes.Therefore, in this minimal configuration, the onset of CSA can be explained by two remaining physical processes: moisture-entrainment-convection (MC) feedback and Convective Heating Overturning Circulation (CHOC) feedback.
The MC feedback involves high detrainment of moist air from deep convective cores.This process locally moistens the surrounding environment, making it more prone to future convection than drier zones (Grabowski & Moncrieff, 2004;Tompkins, 2001).The CHOC feedback, instead, involves deep convection releases latent heat, which amplifies positive buoyancy and pressure perturbations, leading to overturning circulation and convection (D.Yang, 2018).D. Yang (2019) showed that the CHOC feedback could produce CSA even in the absence of the MC feedback, with a much smaller increase in the variance of precipitable water and two different clusters rather than a single moist patch.The relevance of MC feedback in CSA can be measured by its sensitivity to the subgridscale (SGS) turbulence scheme.Indeed, in numerical simulation with coarse horizontal resolution (of the order of 1 km) and explicit convection, mixing processes are largely controlled by the SGS turbulence scheme.At this coarse grid resolution deep convection and turbulence motions inside clouds are not resolved and overturning circulations occur in a laminar mode (Bryan et al., 2003).Moreover at coarse grid resolutions, dissipation processes in cold pools (Grant & van den Heever, 2016) and shallow convection (Janssens et al., 2022), which are fundamental for the development of shallow circulations, may depend considerably on the model's choice.
In this study we address the following questions: • Is mixing and MC feedback necessary for convective self-aggregation to occur?• If so, what are the physical mechanisms by which mixing promotes CSA onset?
• Lastly, can this mechanism account for the relationship between CSA and horizontal grid resolution?
The use of CRMs in the above studies has yielded contrasting results regarding these questions.While no sensitivity to the SGS scheme was found by C. J. Muller and Held (2012) with the SAM (the System for Atmospheric Modeling, M. F. Khairoutdinov & Randall, 2003) model, Tompkins and Semie (2017) found an opposite result with the WRF (the Weather Research and Forecasting Model, Skamarock et al., 2019) model.In particular, they found that a high mixing rate of water vapor (high eddy diffusivity) is a necessary condition for convective organization to occur.This result was also confirmed by the experiments of Shi and Fan (2021).Using Cloud Model 1 (CM1), they demonstrated that different turbulent parametrizations in the PBL can substantially affect the initiation of self-aggregation.
The idealized model of Biagioli and Tompkins (2023) offers one possible solution, where the most important factor for determining the onset of aggregation is found to be the maximum convective-free distance prior to clustering.A significant distance would correspond to a substantial humidity perturbation, which could promote the formation of drier-than-average regions.These regions may expand and lead to CSA, confirming the previous hypothesis by Windmiller and Craig (2019) in that "large enough fluctuation in the humidity content has to be present for self-aggregation to start." To test the ideas derived from those idealized models, this paper thoroughly investigates the impact of mixing and horizontal grid resolution in CSA with different CRMs.The analysis is carried out with the SAM and the WRF models by including the sub-grid scale mixing processes in the Froze Moist Static Energy budget and by studying the horizontal variability of convection through spectral analysis.Energy spectra of quantities such as kinetic energy (KE) and vertical velocity can be used as a measure of the model's ability to reproduce the correct energy statistics (Skamarock, 2004) and to verify the model dynamics against the observations (Lindborg, 1999;Nastrom & Gage, 1985).This is especially true in CSA studies, where the horizontal variability of water vapor and atmospheric motion play a fundamental role.For example, Beucler and Cronin (2019) used a Fourier analysis to see on which spatial scales the different forcing for aggregation occurs and found that longwave radiative fluxes operates at large scales (1,000-5,000 km), while shortwave and surface fluxes operates at smaller scales (500-3,000 km).
Lastly, it is demonstrated here that an efficient updraft dilution, either created by explicit or implicit mixing processes, encourages more large-scale variability of convection (Mapes & Neale, 2011) in the early stages.This process also creates a large enough dry perturbation in the free troposphere to start the radiative feedback loops.We also show how such perturbation varies with different horizontal grid sizes and how it can explain the sensitivity of CSA to horizontal resolution, independent of the amount of low-level cloud.Thus, the representation of small-scale mixing processes becomes fundamental for establishing free-tropospheric (FT) drying, which is considered a prerequisite for the onset of CSA (B.Yang & Tan, 2020;Shamekh et al., 2020;Yanase et al., 2022).
The paper is organized as follows.First, in Section 2, SAM and WRF numerical models are described and the simulations setup, with all the sensitivity experiments, is detailed.The results are divided into four sections.Section 3.1 focuses on the general characteristics of the self-aggregation state and its impact on domain-average statistics.Section 3.2 shows the analysis of power spectra for kinetic energy, humidity and temperature.Section 3.3 studies the relevance of MC feedback in the triggering of CSA.Section 3.4 summarizes the hypothesized mechanism for triggering CSA and compares it with the radiative mechanisms as caused by the WTG hypothesis.The sensitivity of the results to initial conditions, turbulent Pr number and horizontal resolution, is examined in Section 4. The results are discussed in Section 5 and conclusions are given in Section 6.

Numerical Models
The two cloud-resolving models chosen for this work are the System of Atmospheric Modeling (SAM version 6.10.5;M. F. Khairoutdinov & Randall, 2003) and the Weather Research and Forecasting Model (WRF version 4.2.2;Skamarock et al., 2019).Such models offer a solid benchmark for studying the sensitivity of convective aggregation to turbulent and numerical diffusion because:  Muller and Held (2012); Wing and Emanuel (2014); D. Yang (2018); Patrizio and Randall (2019) for SAM); (d) both models participated to the RCEMIP project (Wing et al., 2020).
The physical parametrization (subgrid-scale turbulence, planetary boundary layer, microphysics, radiation, and surface layer) have been fixed to those used by C. J. Muller and Held (2012) for the SAM model and by Tompkins and Semie (2017), (see their smag3dpbl simulation in Table 1) for the WRF model, to make our results directly comparable to their work.All parametrizations used are summarized in Table 1.
Eddy viscosity models do not allow any energy transfer from the SGS scales to the resolved scales.They are purely dissipative and model the turbulent fluxes of momentum (the Reynolds stress term, τ ij = u′ i u′ j , where u′ i is the velocity at subgrid scales) and other scalar quantities, ϕ, in the following form: (1) where the overbar denotes resolved quantities, S ij = 1/2(∂u i /∂x j + ∂u j /∂x i ) is the resolved strain rate tensor and K h,v is the horizontal or vertical eddy viscosity, also called mixing coefficient and Pr is the turbulent Prandtl number.For both SAM and WRF we choose the anisotropic implementation of the Smagorinsky-Lilly closure, where two different mixing coefficients are specified for the horizontal and the vertical direction.
In both SAM and WRF, the definition of K h,v , using the notation by Simon and Chow (2021), takes the following form: where C s is the Smagorinsky constant, l h (l v ) the horizontal (vertical) mixing length scales, and S β is the magnitude of the resolved strain rate tensor modified by a factor based on the local Richardson number.While the ratio between horizontal and vertical mixing coefficient is the same for both models, the main differences between SAM and WRF rely on the calculation of l v , the values of C s and the turbulent Pr number.The latter parameter is a prescribed constant value and it is fundamental not only because it is used in the mixing length calculation (see Text S3 in Supporting Information S1), but also because scalar diffusivity coefficients are simply obtained by dividing the momentum eddy viscosity coefficient by the Pr number.This is valid both for SAM and WRF models.Therefore, at similar values of eddy viscosity, smaller Pr numbers will cause larger scalar diffusivity.(Collins et al., 2006) RRTMG (Iacono et al., 2008) Microphysics SAM1MOM (M.F. Khairoutdinov & Randall, 2003) Purdue Lin (Chen & Sun, 2002) Surface layer Monin-Obukhov similarity Revised MM5 similarity (Jiménez & Dudhia, 2012) Subgrid-scale turbulence 3D Smagorinsky 3D Smagorinsky (Smagorinsky, 1963) PBL None Yonsei University, YSU (Hong et al., 2006) As shown in SAM and WRF large domain simulations are initialized from the equilibrium sounding of a corresponding smaller domain RCE simulation (averaging over the last 20 days).The final RCE equilibrium reached by the small simulations presents some differences between SAM and WRF, even though they start from similar initial profiles.In particular, the WRF profile is moister and warmer than that of SAM and there is a lower Convective Available Potential Energy (CAPE) is present (see Figures S1 and S2 in Supporting Information S1).Similar differences between SAM-CRM and WRF-CRM model were found in RCEMIP project (Wing et al., 2020, see their Figures 7g and 8g).Different initial conditions may strongly influence the sub cloud layer properties and CSA as demonstrated by Cerlini et al. (2023).Further details about smaller domain simulations can be found in Text S3 in Supporting Information S1.
Simulations are divided into three groups: main experiments, sensitivity experiments to initial conditions, sensitivity experiments to Pr number and sensitivity experiments to horizontal resolution.All experiments are listed in Table 2.
The main experiments are performed with a horizontal grid resolution of 3 km.They compare the standard version of the model (standard Smagorinsky and standard numerics) with the version of the model obtained by switching off the Smagorinsky turbulence model.This is obtained by imposing C s = 0 in both SAM and WRF.While in SAM this corresponds to switching off completely the physical mixing (since a minimum coefficient is set on the mixing length prior to the mixing coefficient calculation), in WRF a small constant background mixing is left (see Text S2 in Supporting Information S1).The main experiments are run for at least 100 days and output fields are stored every hour.In the following, the first three letters of the simulation name represent the model name (e.g., SAM, WRF), while the number 0 is appended to simulations with C s = 0 (e.g., SAM0, WRF0).
Updraft entrainment and dilution is influenced not only by the mixing coefficient and the turbulence model but also by the properties of environmental (entrained) air.Therefore we perform an initial conditions sensitivity experiment, by initializing WRF and WRF0 simulations from the same initial conditions of SAM.Such simulation will be named WRFs and WRF0s.This is to prove the robustness of our results with respect to initial conditions.
One way to reduce the numerical mixing inherent to discretization schemes is to use a finer grid resolution.Therefore, the grid resolution sensitivity experiments are designed to investigate the impact of a finer horizontal resolution (1 km) on physical and numerical mixing.Two different high-resolution simulations of WRF are run: WRFh, which is the same as WRF, but with a resolution of 1 km; WRF3h, which is as WRFh, but the Smagorinsky constant is increased by three times.The same experiments are repeated for the SAM model: SAMh, which is the same as SAM, but with a resolution of 1 km; and SAM3h, where the resolution is 1 km, the Smagorinsky constant is increased by a factor of 3 and the Pr is reduced to 1/3 as the WRF value.To obtain a stable solution for the SAM3h experiment, we had to decrease the time step from 10 to 5 s.In the following, simulations with 1 km resolution have a small "h" in the simulation name (WRFh, WRF3h, SAMh, SAM3 h).
High-resolution simulations are run for at least 30 days and output fields are stored every hour.

Convective Aggregation State
Figure 1 shows the 2D snapshots of Outgoing Longwave Radiation (OLR) for the main experiments after 100 days.Convective self-aggregation occurs in WRF, SAM and SAM0: a small circular moist patch with deep convective activity (OLR <120 W m 2 ) of diameter around 200 km is surrounded by a larger and very dry region with suppressed convective activity (OLR >280 W m 2 ).In WRF0, on the other hand, convective activity is homogeneously distributed across the domain without the formation of extremely dry regions (Figure 1c).
The aggregated simulations with active SGS mixing, SAM and WRF, show a similar temporal evolution of domain-average OLR, which increases by almost 50 W m 2 from the beginning of the simulation, reaching a stable value after 40 days (Figure 2).SAM0 shows the same increase but with a slower rate and a stable average OLR is reached after 60 days.WRF0 shows the same initial value of domainaverage OLR throughout the simulation period.
Therefore, the reduction of horizontal diffusion in WRF is found to prevent CSA, as demonstrated by Tompkins and Semie (2017).The same cannot be said for SAM, confirming the previous results of C. J. Muller and Held (2012).SAM0 simulation demonstrates that the removal of turbulent mixing in SAM model, slows down the aggregation process, without preventing it.It is important to note that to characterize the effective slow down of aggregation caused by mixing, we should have run an ensemble of RCE simulations to take account for the stochastic variability of convective selfaggregation.Therefore, our interpretation of the role of mixing will focus on the triggering phase of CSA (the first 10 days) and only on whether a model exhibits or not a final aggregated state, rather than focusing on its degree of organization or its temporal evolution.
Figure 3 shows the initial (0-5 days, black lines) and final (last 10 days, red lines) profiles of temperature and relative humidity in the main experiments.Initially SAM has a lower temperature throughout the troposphere.However, after aggregation, the temperature profile reaches values comparable to WRF.SAM0 is not shown in Figure 3a since as no differences are found with respect to SAM.The increase in temperature in WRF with aggregation is smaller relatively to that of SAM, but it is still evident.WRF0 shows no changes with respect to the initial profile in WRF and is therefore not shown in Figure 3a.In the initial phase, the lower troposphere (below 6 km) is moister in WRF than in SAM and SAM0.No significant difference in relative humidity is found between SAM and SAM0, while the opposite is true for WRF and WRF0.Over the lower troposphere, WRF0 shows a very moist and well mixed profile with respect to WRF.This profile remain stable throughout the simulation period.The upper troposphere (above 6 km) of aggregated simulation (WRF, SAM, SAM0) is drier with respect to WRF0 in the initial stages and shows similar values even comparing different models.After aggregation, the domain dries out and similar profiles of RH are observed in WRF, SAM and SAM0.
The dryness of the upper free troposphere is the only factor common to the aggregated simulations with both models.This factor is already present in the early stages of the simulation as shown in Figure 3b even if the different models starts with different temperature and different lower tropospheric relative humidity profiles.The initial drying of the free troposphere is essential for CSA (Shamekh et al., 2020;B. Yang & Tan, 2020;Yanase et al., 2022), since it allows a more efficient radiative cooling of the boundary layer in the dry regions and the onset of a shallow circulation that transports MSE from dry to moist regions (upgradient), contributing positively to CSA.This drying has to be caused by a relatively fast and efficient process which is acting on the triggering phase of CSA, which is considered here to occur during the first 5 days.In order to better characterize this dry perturbation we will analyze the horizontal variability of convection through the energy spectra (Section 3.2) and then focus on the physical processes responsible for creating such a perturbation (Section 3.3).
We obtain contrasting results between SAM and WRF models by switching off the SGS turbulence model.We hypothesize that this difference rely on the different implicit numerical mixing between the two models.This hypothesis will be demonstrated later by the sensitivity experiments performed on finer horizontal resolution.

Energy Spectra and Horizontal Variability of Convection
We compute energy spectra for horizontal Kinetic Energy (KE), vertical velocity (w), and virtual potential temperature perturbation (θ′ v where θ v = θ (1 + 0.608q v )) both in the free troposphere (averaging from 3 to 10 km) and in the boundary layer (averaging up to 2 km).Then we also compute energy spectra for the Total Water Path (TWP, calculated as the sum of all water species), integrating either over the whole column or over the boundary layer  (from 0 to 2 km).While KE and w variables are important to understand the horizontal variability of atmospheric motions, θ′ v and TWP give us a measure of buoyancy and humidity fluctuations.The results for the free troposphere (and the whole column TWP) are shown in Figure 4, while the boundary layer calculation is reported in Figure S3 in Supporting Information S1.The horizontal one-dimensional PSD of each variable is obtained following the procedure of D. Durran et al. (2017).The PSD is then multiplied by the corresponding wavenumber in order to have a direct correspondence between the variance of the variables and the areas underneath the reported spectrum curves.All the spectrum figures contain two reference lines (dashed gray lines): the vertical line indicates wavelengths equal to 12 km which is taken as an indication of the effective resolution of the model (4Δx, D. R. Durran, 2010); the oblique line shows the k 5/3 power law which is commonly observed where there is a turbulent energy cascade.Lai and Waite (2023) found that RCE simulations with the WRF model at 4 km horizontal resolution can show the 5/3 power low in the KE spectrum in the upper troposphere.However, from their spectral budget analysis, they found that this spectrum to be generated by a balance between forcing by buoyancy flux and removal by vertical energy flux, instead of a turbulent cascade (where nonlinear advective transfer is the main contribution to the budget analysis).
There are two main differences between the kinetic energy distribution of SAM and WRF as shown in Figures 4a  and 4b: (a) SAM atmospheric motions, both in the vertical and in the horizontal, contain more energy across all .Power spectral densities of (a) horizontal Kinetic Energy (KE), (b) vertical velocity (w), (c) perturbation virtual potential temperature (θ′ v ), averaged over the free troposphere (between 3 and 10 km) and (d) total water path (TWP) averaged over the whole column.The reported values are time-averaged in the initial 5 days of simulation (black lines) and the last 10 days (red lines).The PSD is multiplied by the corresponding wavenumber to have a direct correspondence between the variance of the variables and the areas underneath the curves.The oblique gray dashed line represents the k 5/3 power law, while the vertical gray dashed line marks the effective resolution of the model, taken as 4Δx.
scales; (b) the difference between kinetic energy spectra increases at smaller scales where WRF is found to damp small fluctuations (wavelengths less than 4Δx), whereas in SAM small horizontal and vertical motions shows a large small-scale variance in the total kinetic energy.These two diverging behavior of kinetic energy reflect different modeling approaches at short wavelengths: WRF uses explicit filters to add scale-selective dissipation and damp the shortest wavelengths which are not well resolved by the model grid, while SAM does not apply any explicit scale-selective filter allowing a larger energy input at the smallest resolvable scales.
One interesting aspect is that the decrease of turbulent mixing in SAM0 did not significantly alter the kinetic energy spectrum both in the horizontal and in the vertical.This may be an indication that in SAM, at grid resolution larger than 1 km, the energy dissipation is controlled mainly by the numerical discretization.
The development of CSA in WRF and SAM simulations causes similar changes to the distribution of energy in both WRF and SAM.With CSA, the spectrum of horizontal KE and vertical velocity becomes more dominated by large-scale motions rather than small-scale motions as found in the non-aggregated case (Figures 4a and 4b).Such an increase of large-scale variance is not observed in the WRF0 simulation, whose spectrum remains constant throughout the simulation.In particular, WRF0 is characterized by smaller and less energetic horizontal and vertical motions already from the beginning of the simulation.Indeed an increase of energy at the smallest resolvable scales is evident in 4a.
The diverging behavior at small scales between SAM and WRF is also visible in the spectra of TWP and θ′ v as shown in Figures 4c and 4d.Larger humidity and buoyancy fluctuations are present in SAM at small scales.This was evident also by the large presence of small cloud structures in SAM simulations as shown by the OLR snapshots (e.g., Figures 1a and 1b).The energy peak of TWP for WRF is found at wavelengths around 20 km.For SAM, this peak is lower and it is found at slightly larger wavelengths (about 30 km).However, the TWP spectrum is quite flat from wavelength smaller than 30 km, indicating that strong updrafts and cloud structures in SAM can be found at very different scales.Instead, in WRF, as shown also in Figure 1a, small deep convective cores and clouds are rare.
When CSA occurs (red lines in Figures 4c and 4d), there is an increase in humidity and virtual potential temperature variance at large scales.Such increase is very large for TWP.The same production of mesoscale humidity variance was observed by Yanase et al. (2022) and Janssens et al. (2022) in the aggregation of deep convection in several models and by C. Bretherton and Blossey (2017) in the mesoscale aggregation of shallow cumulus (see their Figure 16).The most prominent change of SAM and WRF with the reduction of turbulent mixing is evident in the initial stages in the variability of TWP, as shown in Figure 4d.Lower horizontal turbulent mixing in SAM0 and WRF0 causes a smaller variance of TWP across all scales.In this way, smaller perturbations of humidity will delay the onset of aggregation (as in SAM0), or prevent it (as in WRF0).The reduction of humidity variance in WRF is greater than in SAM.
The reduction in free tropospheric moisture variance with decreased horizontal mixing is well correlated with a more active boundary layer in both SAM and WRF as seen in the vertical velocity spectrum in Figure S3b in Supporting Information S1.

FMSE Variance Budget
The time evolution of CSA is characterized by an increase in the Frozen Moist Static Energy (FMSE) variance, where FMSE is defined as: where c P is the specific heat of dry air, T is the air temperature, g is gravitational acceleration, z is the height above the surface, L v is the latent heat of vapourization, q v is the specific humidity with respect to water vapor, L f is the latent heat of fusion, and q i is the specific humidity with respect to ice condensates.The evolution of its horizontal spatial variance, ĥ′2 , is shown in Figure 5a for the main experiments.(For each quantity x, we denote "x" as its density-weighted vertical integral ∫ z top 0 xρdz.The horizontal domain mean of x is instead denoted as {x} and the anomaly as x'.).WRF and SAM exhibit a similar exponential growth of variance and they reach a stable equilibrium state after 40 days.The final variance values are also very similar between the two models.In order to study the feedback that cause the increase of the FMSE variance, we perform the budget analysis according to the evolution equation of ĥ′2 by Wing and Emanuel (2014): where SEF are the total surface enthalpy fluxes (the sum of the latent and the sensible heat flux), NetSW and NetLW are the column shortwave and longwave radiative flux convergence, and ∇ h ⋅ û → h is the horizontal divergence of the vertically integrated flux of h, which is evaluated as a residual of the remaining terms.Radiative processes have very similar magnitude and evolution between SAM and WRF: they contribute positively to self-aggregation throughout the simulation time.In particular, longwave feedback are the drivers of CSA in the first stages up to 10 days, with WRF having higher correlations than SAM.After 10-20 days, longwave and shortwave feedback have similar magnitudes and remain almost constant, indicating that they are responsible for maintaining CSA.Surface enthalpy flux contributions are positive up to 20 days for both SAM and WRF.However, in SAM they exert stronger positive feedback than in WRF, especially in the early stages where their magnitude is comparable with that of longwave radiative processes.After 20 days, the contribution of surface fluxes becomes negative, thereby opposing self-aggregation.Then, in SAM they return to being slightly positive, while in WRF they remain slightly negative.This is one of the few differences between WRF and SAM simulations, regarding this budget analysis.Advective contributions are counter to CSA in the early stages, but they become slightly positive in a time window of 10 days, which corresponds to the time when the convection starts to organize into a single cluster.This time window starts earlier in SAM than in WRF and overlaps with the time when the surface fluxes feedback becomes slightly negative.After the convective cluster has formed, the advective processes return to be negative, counteracting the CSA and the radiative feedback that sustain it.
From the above FMSE analysis, convective self-aggregation is primarily triggered by radiative and surface fluxes feedback.However the contribution of SGS turbulent mixing is not taken into account in Equation 6.By using Equation 2 for modeling turbulent diffusion of FMSE, one can add the SGS mixing contribution to the FMSE variance budget as follow: Figure 6 offer a zoom of Figure 5 during the first 10 days, including also the effect of SGS mixing.The increase in the FMSE variance in Figure 6a is quite similar between WRF and SAM, while it is slower in SAM0 and absent in WRF0.This correlates with the total diabatic feedback, as they show similar increase in magnitude for SAM and WRF, while a slower increase in SAM0 and no increase at all for WRF0.Surprisingly, the SGS mixing term shows much larger positive correlations than the total diabatic feedback in the very early stages for SAM and WRF, implying that it strongly favors CSA and it cannot be neglected at such horizontal resolution.Looking only at SAM and WRF simulation, this result implies that a strong cooperation between SGS mixing and diabatic feedback is necessary to start the diabatic feedback loop which is then responsible for the expansion and maintenance of dry patches.
However, switching off such feedback in SAM0 does not prevent a constant increase in the magnitude of diabatic feedback, while it does in WRF0 where the amplitude of diabatic feedback remain constant.We hypothesize that this behavior is related to the implicit numerical entrainment which is present in SAM, when used at coarse resolutions.We will test this hypothesis by looking at the FMSE budget at finer resolution with different explicit and implicit mixing.

Triggering Mechanism of CSA
From the analysis of the horizontal variability of convection and the FMSE budget, it is clear that mixing processes are fundamental for triggering CSA.In particular, we argue that the radiative feedback loop which is responsible for the initial formation and the expansion of dry patches in the very initial stages, is initiated through updraft dilution in a dry environment instead of the radiative subsidence due to the WTG.This is consistent with previous findings by Yanase et al. (2022).Following the work by Tompkins and Semie (2017), we suggest that the free-tropospheric drying and the associated radiative subsidence are strictly dependent on lateral mixing and updraft dilution.Models with an efficient mixing, either numerical or explicit, will allow for greater large-scale variability of convection and therefore greater dry perturbations over the free troposphere.Only when a sufficiently strong dry perturbation is established by updraft dilution and reduced convective heating and moistening, then a strong boundary layer radiative cooling in dry region is able to start the upgradient MSE shallow circulation and expand dry patches, as also found by Shamekh et al. (2020).Figure 7 shows respectively radiative cooling, static stability, cloud fraction, radiative velocity (calculated as the ratio between radiative cooling and static stability following Bony et al. ( 2016)) and the actual velocity in the driest regions (the lowest quartile of blockaveraged Column Relative Humidity, CRH) averaged over the first 5 days.The strongest differences in the radiative cooling profiles are evident in the boundary layer below 2 km (see Figure 7a).Simulations with convective self-aggregation shows an average cooling of 1.5 K/day over the boundary layer with a peak at 1 km corresponding to the maximum low-level cloud amount (see Figure 7c).However, this cooling is not present in WRF0.The reason behind the absence of such cooling cannot be traced back to the decrease of radiative subsidence, since Figure 7d shows that between 2 and 10 km the radiative velocity between WRF and WRF0 is quite similar.Instead this difference has to be related to the high amount of the mid-level clouds (Figure 7c) and the positive low-level actual velocity (Figure 7e) in the driest region, which implies that convection is able to penetrate into the free troposphere of dry regions and destroy any nascent dry perturbation.
Although the WTG path, radiative velocity and clear-sky convergence may be useful for explaining different anvil cloud fraction, cloud top and temperature, Figure 7 shows that they cannot explain the different sensitivity of CSA in SAM and WRF to mixing processes.In fact, despite the different radiative cooling in the upper troposphere and the different static stability (Figures 7a and 7b), the WTG velocity, as diagnosed from the radiative cooling, is quite similar between SAM and WRF between 2 and 8 km.This is because the warmer atmosphere in WRF cause a larger radiative cooling which compensates for the greater static stability with respect to SAM.A similar effect was found also by Shamekh et al. (2020) by varying SST.

Sensitivity to Initial Conditions
Entrainment and updraft dilution depend also on the ambient temperature and relative humidity.A drier atmosphere, especially in the lower troposphere such as that of SAM, would amplify the dilution effect.On the other hand, in warmer atmospheres, such as that of WRF, entrainment is more effective in reducing updraft buoyancy (Singh & O'Gorman, 2013).Different mixing efficiency due to different environment could lead to different thermal stratification (static stability).Through the Weak Temperature Gradient hypothesis this would directly impact the radiatively driven subsidence (and all the deep convective circulation in general), the clear-sky convergence and also the anvil cloud fraction (Bony et al., 2016) with possible consequences for convective self-aggregation.However, based on our previous results, we have shown that the free-tropospheric drying necessary to start the radiative feedback loop, is not directly caused by the WTG path, but by the initial dry perturbations set by mixing processes.Therefore we do not expect self-aggregation in WRF0 even if we start the model with a colder atmosphere and a drier lower troposphere.
Figure 8 further confirm our findings.Even by starting WRF and WRF0 by using SAM initial RCE profiles (respectively WRFs and WRF0s), the relative occurrence of convective-self aggregation in the two simulations is not changed.Indeed after 10 days, Figure 9 shows that the temperature and humidity profiles of WRFs and WRF0s are very close to their corresponding simulation initialized with warmer and moister atmospheres (WRF and WRF0).A small free tropospheric drying is evident in WRF0s (see Figure 9b), but this is not sufficient to start the radiative feedback loop.
The WRF model seems to be very efficient in moistening the lower troposphere with respect to SAM, even with a very small amount of mixing, as in WRF0s.This is mainly due to the convectively induced moistening which is favored by a reduced updraft dilution, as denoted by the difference between WRFs and WRF0s in Figure 9b.Moreover, it could also be related to the fact that SAM does not have a planetary boundary layer scheme, while WRF does.However, the impact of the planetary boundary layer scheme in WRF has already investigated by Tompkins and Semie (2017), and they found that including vertical mixing with a PBL scheme favors CSA.Therefore the results obtained for WRF0s should be conservative with respect to the choice of the PBL scheme.
Figure 10 shows the correspondent increase in diabatic feedback and SGS mixing feedback for WRFs with respect to WRF.Another difference is the  strong peak observed in WRFs in the FMSE variance at 12 hr with respect to that of WRF.This peak is the consequence of a larger effective updraft dilution in dry regions in WRFs with respect to WRF.The reduction in the number of deep convective cores is reported in Figure S4 in Supporting Information S1.
We conclude that the free-tropospheric drying caused by updraft dilution in WRF depends mainly on the SGS mixing and that a drier environment is not sufficient to allow convective self-aggregation to develop in the absence of lateral mixing.This necessary condition for CSA triggering is again shown in Figure 11a, where we highlight the difference in large-scale horizontal variability of convection between WRFs and WRF0s during the first 5 days.
Figure 12a shows the resulting differences in cloud fraction profiles.The most significant difference between WRF and WRFs is an increase in the height at which low-level cloud form and a decrease of the height of anvil clouds.WRF0s and WRFs generally show a larger anvil but they have same distribution of WRF and WRF0 in the lower troposphere.The decrease in anvil cloud height and their increase in coverage can be traced back to the colder atmosphere of WRFs and WRF0s in agreement with the study of Stauffer and Wing (2022), and it also partly explain the difference between SAM and WRF simulations.However, the large differences in upper tropospheric and low-level cloud cover profile between SAM and WRFs depend mainly on the microphysics scheme.

Sensitivity to Pr Number
By default, the WRF and SAM models use different Pr numbers, equal to 1/3 for WRF and 1 for SAM.The Pr number is a fundamental constant in the calculation of the eddy viscosity and eddy diffusivity.First, it appears in the calculation of the resolved strain tensor (see Text S2 in Supporting Information S1) as a weighting factor in the reduction of mixing in very stable atmospheric conditions.Therefore, holding the resolved strain tensor and the buoyancy frequency constant (hence a constant Richardson number), the smaller the Pr number, the greater the reduction of mixing in a stable stratified state.Secondly, the eddy diffusivity is derived by dividing the eddy viscosity by the turbulent Pr number.Therefore a larger Pr number implies a smaller eddy diffusivity, at constant eddy viscosity.
Due to its influence in scalar mixing processes, we expect a large sensitivity of CSA to the Pr number.In particular, given Equation 6, the smaller the Pr number, the larger the contribution of SGS mixing in the increase of FMSE variance and the larger the moisture-convection feedback that allows free tropospheric drying.
Figure 13 shows the OLR field after 20 days.The onset of CSA is clearly visible in SAMPr03, while WRFPr1 shows no sign of CSA.Therefore, increasing the Pr number to 1 in WRF is enough to prevent CSA, as also found by Shi and Fan (2021) for the CM1 model.After 100 days, WRFPr1 exhibits always random convection like WRF0 and WRF0s, while SAMPr03 aggregated state is very similar to SAM and SAM0 simulations.
As expected, decreasing the Pr number in SAM increases the strength of SGS mixing and diabatic feedback as shown in Figures 14b and 14c.This corresponds to a faster increase in the FMSE variance (Figure 14a).
The variability of the Total Water Path in SAMPr03 decreases at all scales (Figure 11b).However the large-scale variability remains always bigger than WRF0 and of the order of magnitude of that of SAM and WRF.The same cannot be said for WRFPr1 which instead exhibits a smaller large scale variability very similar to that of WRF0.The similarity between WRF0 and WRFPr1, SAM and SAMPr03 can be seen also by looking at the relative humidity vertical profiles shown in Figure S6 in Supporting Information S1.
Decreasing the Pr number in SAMPr03 doubles the average low-level cloud fraction, reaching values and depths more similar to that of WRF (see Figure 12b).The same increase is found also in the driest regions (not shown).

Sensitivity to Horizontal Resolution
In order to demonstrate that the SAM model has a large numerical mixing that makes CSA insensitive to SGS mixing at coarse resolution of 3 km, we run the main experiments WRF and SAM at finer resolution.We then reintroduce mixing explicitly, through the SGS turbulence model, in order to obtain CSA with both models even at higher resolution.This demonstrates that the total amount of mixing, either implicit or explicit, regulates the large scale variability of convection and therefore the free-troposphere drying necessary to trigger CSA.A correct rescaling of mixing coefficient with the horizontal resolution is needed if we want to maintain a constant largescale variability.
Figure 15 shows the OLR snapshots for finer resolution simulations after 30 days.WRFh and SAMh simulations show no sign of CSA, although a larger scale variability of convection is evident in SAMh compared to WRFh (Figures 15a and 15b).Therefore, reducing the mixing length by a factor 1/3 and reducing the numerical mixing is enough to prevent CSA for both models.CSA is restored by reintroducing explicit mixing in WRF3h (simply increasing by a factor 3 the Smagorinsky constant) and in SAM3h (by decreasing Pr number by a factor 3 and increasing the Smagorinsky constant by a factor 3, to have similar constant to WRF3h).
WRF3h exhibits the strongest increase in FMSE variance, while SAM3h exhibits similar increase to those of SAM and WRF (Figure 16a).The effect of SGS mixing is larger in WRF3h compared to that of SAM3h, which is quite similar to that of all other simulations, except for SAMh, which has a very low correlation for this term.This shows the contributions to the FMSE budget by SGS mixing, diabatic feedback (longwave, shortwave and surface flux) and advection (evaluated as a residual term), respectively.Each contribution has an hourly time-step and a 6-hr running average is applied to the horizontal convergence term.
is another signal of the larger contribution of numerical mixing in the SAM model.Indeed, at the same resolution and with same Pr numbers, WRF needs a much larger explicit mixing to trigger CSA.
In the literature, the dependence of the occurrence of CSA with the horizontal grid resolution is commonly associated to the low-level cloud amount: the finer the resolution, the smaller the low-level cloud amount and the associated low-level radiative cooling which is necessary to start the shallow MSE upgradient circulation (C.J. Muller & Held, 2012;Yanase et al., 2020).Figure 17a shows that this dependence is found for the WRF model, since WRF and WRF3h have larger amount of low-level cloud fraction than WRFh.However, we recall that the main difference between WRF and WRF0 was not in the low-level cloud amount, but on the mid-and upper-level cloud amount in the driest regions (see Figure 7c).On the other hand, Figure 17a shows that the two simulations SAMh and SAM3h have similar low-level cloud fractions on average, demonstrating that we can obtain CSA even with a small amount of low-level clouds.From this analysis we conclude that the average amount of lowlevel clouds does not affect the triggering phase of CSA, while it may be fundamental for its intensification and development.
The most important factor that binds all the aggregated simulations together in the very early stages of the simulation is the large-scale variability of the convection.This is shown in Figure 17b, where the difference in large-scale variance of TWP between the non-aggregated and aggregated cases is highlighted.It is interesting to note the larger large-scale variance of SAMh with respect to WRFh.This difference may be associated to the larger numerical mixing of SAM model with respect to the WRF model.This also reflects the differences found in the horizontal variability of OLR between WRFh and SAMh (Figures 15a and 15b): WRFh has a more random and uniform convection, with smaller dry areas, while SAMh exhibits larger clusters of deep convective cores with larger dry areas.
The spectral analysis results in Figure 17 indicate the importance of the average distance between deep convective cores and the initial size of humidity perturbations, as discovered by the idealized model of Biagioli and Tompkins (2023).Increasing the Smagorinsky constant enlarges the area of influence of each deep convective core, as well as their average size.This increases the mean distance between different cores, resulting in a larger amplitude of dry perturbation introduced in the free troposphere.By maintaining the original updraft width, even at a finer grid resolution, such modification restores CSA.Figures 18a and 18b show the average and maximum convective distance between updrafts over the first 5 days of simulations.Updrafts are detected by selecting the grid cells with a column relative humidity exceeding 0.9.The distance free from convection is calculated using the method developed by Biagioli and Tompkins (2023).The distance values indicate a clear separation between simulations that exhibit CSA and those where convection remains disaggregated.The average and maximum distance free from convection in the aggregated simulations is three times larger than in simulations without CSA.This is evident even just a few hours after the simulation starts (refer to Figure S6 in Supporting Information S1).
This analysis demonstrates a strong correlation between the large-scale variability of convection as measured from horizontal spectra of TWP and that measured from the convective distance metric.The results indicate that the larger the horizontal variability of TWP at large scales, the larger the free-convective distance.Therefore, these sensitivity experiments confirm the results of the idealized model from Biagioli and Tompkins (2023) and the hypothesis formulated by Windmiller and Craig (2019) that self-aggregation requires a significant fluctuation in humidity content to start.It is concluded that coarse horizontal grid resolutions increase this initial fluctuation, either through larger SGS mixing or numerical mixing.

Discussion
This study shows the relevant role of mixing and entrainment on the self-aggregation phenomenon.The level of mixing and entrainment reproduced by numerical simulations is provided by the explicit turbulence model and the numerical diffusion.Accordingly, we have shown a strong sensitivity of CSA to the different turbulence models and spatial resolutions adopted in SAM and WRF atmospheric models.We try to capture the essential aspects here.
On one hand, in the SAM model, a lower turbulent mixing coefficient, hence a weaker MC feedback, does not affect the onset of CSA.Therefore, looking only at this model, modifying the sensitivity of convection to water vapor, cannot by itself destabilize the RCE, as suggested by Emanuel et al. (2014).On the other hand, in the WRF model, turbulent mixing appears to be a necessary condition for CSA to occur.This result is consistent with previous work by Tompkins and Semie (2017).
These two contrasting results can be reconciled by considering all scalar mixing processes either explicit (turbulence, hyperviscosity, or other explicit numerical filters) or implicit (diffusion generated by the advection schemes).Implicit numerical diffusion can provide a large contribution to the total mixing, especially at coarse grid resolution (such as those used in our simulations), and can cause numerical entrainment (Yamaguchi et al., 2011).Therefore, if the SAM model at coarse resolution is characterized by a large numerical diffusion, its solutions will be substantially unaltered by switching off the explicit turbulence model.This approach resembles the one used in Implicit Large Eddy Simulations (ILES).In contrast, in the WRF model, which is built with higher-order numerics for running at coarse resolution, most part of the mixing is provided by the turbulent subgrid-scale parametrization.Accordingly, a substantial change of the WRF solution is observed by switching off the turbulence model.
We have shown that reducing the numerical and turbulent diffusion in SAMh (by refining horizontal resolution) is sufficient to prevent the onset of CSA.A finer resolution and a lower diffusion is found to decrease also the amount of shallow clouds in the driest regions and this may be interpreted as the main cause behind the absence of CSA at finer resolutions (C.J. Muller & Held, 2012).However, the reintroduction of a large explicit mixing in SAM3h allows the onset of CSA at finer resolution, even with a small amount of shallow clouds.Therefore, shallow clouds and their effect on radiative cooling, do not seem to play a significant role in the triggering of CSA.Instead, they may be considered fundamental at later stages for its development and intensification.
The key ingredient common to all numerical simulations showing the onset of CSA is the initial large-scale horizontal variability of convection.This is examined by looking at the horizontal energy spectra of TWP.It is shown that a larger horizontal variability at large-scales is necessary to introduce large dry perturbations in the free troposphere and to cause the initial free-tropospheric drying.Such an initial large-scale variability is found from the first days of the simulation.Therefore, it cannot be created through a slow process such as the radiative subsidence, as calculated by the WTG hypothesis, but must be traced back to a fast and efficient process linked to the environment for convection, such as lateral mixing.An efficient lateral mixing will dilute updraft in dry regions, reducing the convectively induced moistening of the free-troposphere and enhancing the large-scale variability of convection.Once such pre-condition of free tropospheric drying is verified, then a strong deep circulation between moist and dry regions can develop, intruding dry air into the boundary layer and creating a high surface pressure anomaly which is recognized to be fundamental for the development of CSA (Shamekh et al., 2020;B. Yang & Tan, 2020).
The necessity of an initial large scale variability in moisture and convection to trigger CSA was hypothesized by Windmiller and Craig (2019) and it is consistent with recent results by Yanase et al. (2022) and Biagioli and Tompkins (2023).In particular, Yanase et al. (2022) demonstrated with the SCALE model that the radiative mechanism is not caused directly by the WTG velocity, but by the indirect path through the environment for convection.The sensitivity of CSA to diffusion was attributed to the effect of humidity entrainment into updrafts by Tompkins and Semie (2017).Biagioli and Tompkins (2023) offer a reinterpretation of their results by attributing it to the mean updraft size.Here we demonstrate that in CRMs, the initial large-scale variability of convection is regulated by small-scale mixing processes.The importance of energy dissipation in convective permitting models at coarse resolution was also found to be critical for shallow convective self-aggregation by Janssens et al. (2022).Interestingly, the appearance of larger and more energetic structures by damping small-  1).A running average with a 24 hr window is used to smooth out the diurnal cycle.
scale structures (imposing larger values of the Smagorinsky constant) is also a common feature in Large Eddy Simulations of turbulent channel flows (Hwang & Cossu, 2010).
In this paper, the MC feedback was found to be fundamental for the triggering of CSA.Our results appear to contrast with those of D. Yang (2019), who show that CSA can be observed even by removing the MC feedback (homogenizing the water vapor field each 3 hr).However, D. Yang (2019) uses the SAM model at 2 km horizontal resolution, which we have shown to be strongly affected by numerical mixing processes.Moreover, the homogenization time scale (3 hr) of the water vapor could be too large to effectively turn off the MC feedback, as also estimated by Ramirez Reyes (2023).Therefore the final aggregated state obtained by D. Yang (2019), which exhibit a very small variance of precipitable water with respect to what is usually observed in CSA studies, could be model dependent.Further work is needed to replicate the experiments of D. Yang (2019) in other CRMs and to find out which part of the SAM numerics is providing the largest contribution to numerical diffusion when using such model at coarse resolution.This may be the large integration time step in the AB3 scheme, the SOC momentum advection scheme (see Text S1 in Supporting Information S1), or the scalar advection scheme (although this work adopted the highest-order available, Yamaguchi et al., 2011).
This pre-condition on initial large-scale variability of convection for the triggering of CSA, can also explain three different properties that characterize self-aggregation in numerical simulations: (a) the sensitivity to horizontal grid resolution; (b) the sensitivity to domain size, (c) the hysteresis of self-aggregation.Coarse resolutions impose larger mixing (SGS or numerical), which makes updrafts more sensitive to the dryness of the free troposphere, increasing the large scale variability of convection and further increasing the free tropospheric drying.Therefore, at finer resolutions, unless we do not keep fixed the mean updraft size and entrainment, CSA will not occur.At the same time, the larger the domain size, the larger the maximum allowed horizontal variability of water vapor and the larger the possibility to trigger CSA, as also found by Yanase et al. (2022).These two arguments can explain why the regime diagram of CSA (the sensitivity of CSA to horizontal resolution and domain size) is different in different models.Using the SCALE model, Yanase et al. (2020) found CSA to occur even at 500 m, while for the SAM model no self-aggregation was found below 2 km (C.J. Muller & Held, 2012).Regarding the hysteresis of self-aggregation, starting the simulation from an already aggregated condition, sets the initial large-scale variability to a fixed large value and we expect no sensitivity to horizontal resolution (C.J. Muller & Held, 2012) or turbulent mixing (Tompkins & Semie, 2017).
Many studies investigate the role of other important factors in the triggering of CSA besides the subgrid-scale mixing, such as vertical resolution (Jenney et al., 2023), grid anisotropy (De Roode et al., 2022), microphysics (Shi & Fan, 2021), rotation (Carstens & Wing, 2022) and SST (Coppin & Bony, 2015;M. Khairoutdinov & Emanuel, 2013;Wing & Emanuel, 2014).Future work is needed to investigate whether the hypothesized precondition on large-scale horizontal variability in convection in the triggering of CSA would be valid for different values of the above cited factors.

Conclusions
The representation of mixing at small scales remains a substantial problem for many large-scale geophysical flows (Mapes & Neale, 2011).In particular, mixing processes affect the amount of cloud condensate (Jeevanjee & Zhou, 2022), the free-tropospheric relative humidity in the tropics (Grabowski & Moncrieff, 2004), and the organization of convection (Takemi & Rotunno, 2003;Tompkins & Semie, 2017).Therefore small-scale mixing processes may have important implications for the climate.
This work focuses on convective self-aggregation, which occurs spontaneously in idealized simulations of RCE.The occurrence of this phenomenon is sensitive to mixing and entrainment processes, as represented by the explicit turbulence model and the numerical diffusion.During the initial stages, prior to the onset of CSA, increased mixing makes updrafts more sensitive to the dryness of the free troposphere, which strengthens the moisture-convection feedback.Efficient lateral mixing dilutes updrafts in dry regions and enlarges them in moist regions, resulting in an increased distance between convective cores and a stronger humidity perturbation in the free troposphere.The study finds that a significant variability of free tropospheric moisture is required for the development of the radiative feedback loop and the occurrence of CSA, which supports the previous findings and hypotheses of Biagioli and Tompkins (2023) and Windmiller and Craig (2019) using idealized models.When refining the horizontal grid resolution and decreasing numerical diffusion, it is necessary to ensure enough turbulent mixing to maintain a large variability of free tropospheric moisture and allow for the spontaneous development of convective aggregation.In these conditions, convective self-aggregation can develop even with a relatively small amount of shallow clouds, indicating that they are not a key factor in determining the dependence of CSA on horizontal resolution.
Until horizontal grid resolution down to the inertial range for deep convection is achieved (100 m, Bryan et al., 2003), non-homogeneous entrainment processes will not be resolved and the response of deep convection to moisture perturbation will be tightly linked to the numerics and the SGS turbulence parametrization adopted by the models.Less idealized Global Cloud Resolving Models (GCRM) should take into account such dependence of mesoscale organization both on numerical and SGS mixing.TWP energy spectra have been shown to be a useful diagnostic tool for assessing the initial large scale variability of convection and its influence on mesoscale organization.In particular, such large-scale variability should not vary with the coarsening/refinement of horizontal resolution and should be comparable to observations.A correct representation of large-scale horizontal variability of convection will also affect the reproduction of MJO events, as demonstrated by Holloway et al. (2013).
In pursuit of these results, new families of turbulence models should be adopted for anisotropic and coarse resolution grids (Cimarelli et al., 2019;Honnert et al., 2021) and the idealized RCE simulations could provide a useful setting for studying their impact on deep convection.
(a) they are widely used across the scientific community; (b) they have different mathematical formulations, different numerics and different physics; (c) the mechanisms behind convective self-aggregation have already been investigated for both models (e.g., see Tompkins and Semie (2017); Colin et al. (2019); B. Yang and Tan (2020), for WRF and C. S. Bretherton et al. (2005); C. J.

Figure 2 .
Figure 2. Domain averaged OLR time evolution for the main experiments.

Figure 3 .
Figure3.Vertical profiles of (a) absolute temperature and (b) relative humidity for the main experiments.Black lines show the average on the first 5 days, while red lines show the average over the last 10 days.WRF0 and SAM0 temperature profiles are not reported in (a), since in the initial stage they are identical to WRF and SAM profiles, respectively.In the final stages, WRF0 does not show any change, while SAM0 has the same profile as SAM.

Figure 4
Figure 4. Power spectral densities of (a) horizontal Kinetic Energy (KE), (b) vertical velocity (w), (c) perturbation virtual potential temperature (θ′ v ), averaged over the free troposphere (between 3 and 10 km) and (d) total water path (TWP) averaged over the whole column.The reported values are time-averaged in the initial 5 days of simulation (black lines) and the last 10 days (red lines).The PSD is multiplied by the corresponding wavenumber to have a direct correspondence between the variance of the variables and the areas underneath the curves.The oblique gray dashed line represents the k 5/3 power law, while the vertical gray dashed line marks the effective resolution of the model, taken as 4Δx.

Figure 5 .
Figure 5.The budget of FMSE variance for all 100 simulations including: (a) time evolution of the domain mean spatial variance of the vertical-integrated Frozen Moist Static Energy (FMSE) (J 2 /m 4 ) for the main experiments; (b-e) shows the contributions to the FMSE budget by longwave, shortwave, surface flux and advection (evaluated as a residual term), respectively.Each contribution is averaged daily and a 5-day running average is applied to the horizontal convergence term.

Figure 6 .
Figure6.The budget of FMSE variance for the first 10 days including: (a) time evolution of the domain mean spatial variance of the vertical-integrated Frozen Moist Static Energy (FMSE) (J 2 /m 4 ) for the main experiments; (b-e) shows the contributions to the FMSE budget by SGS mixing, diabatic feedback (longwave, shortwave and surface flux) and advection (evaluated as a residual term), respectively.Each contribution has an hourly time-step and a 6-hr running average is applied to the horizontal convergence term.

Figure 7 .
Figure 7. Vertical profiles of (a) radiative cooling, (b) static stability, (c) cloud fraction, (d) radiative velocity (obtained as the ratio between radiative cooling and static stability), (e) actual velocity.All quantities are averaged over the first 5 days and over the driest regions.Following the approach by C. S.Bretherton et al. (2005), all quantities are sorted into four quartiles by the block-averaged Column Relative Humidity (CRH), computed by dividing the simulations domain into blocks with area of 96 km 2 .Driest regions correspond to the lowest quartile.

Figure 9 .
Figure 9. Vertical profiles of (a) absolute temperature and (b) relative humidity for the simulation WRF, SAM, WRFs and WRF0s.Black lines show the average on the first 6 hr in order to show the corresponding initial conditions between SAM, WRFs and WRF0s and their difference with respect to WRF.Red lines show the profiles averaged over day 10.After 10 days the profile of temperature and relative humidity of WRFs and WRF0s are very close to that of the respective simulations WRF and WRF0, indicating a small influence of initial conditions on the occurrence of selfaggregation.

Figure 10 .
Figure10.The budget of FMSE variance for the first 5 days including: (a) time evolution of the domain mean spatial variance of the vertical-integrated Frozen Moist Static Energy (FMSE) (J 2 /m 4 ) for the initial sensitivity experiments; (b-e) shows the contributions to the FMSE budget by SGS mixing, diabatic feedback (longwave, shortwave and surface flux) and advection (evaluated as a residual term), respectively.Each contribution has an hourly time-step and a 6-hr running average is applied to the horizontal convergence term.

Figure 11 .
Figure11.Power spectral densities of total water path (TWP) averaged over the whole column and over the first 5 days for (a) sensitivity experiments on initial conditions, (b) sensitivity experiments to Pr number.The PSD is multiplied by the corresponding wavenumber to have a direct correspondence between the variance of the variables and the areas underneath the curves.The oblique gray dashed line represents the k 5/3 power law, while the vertical gray dashed line marks the effective resolution of the model, taken as 4Δx.

Figure 12 .
Figure 12.Vertical profiles of domain-mean cloud fraction avaraged over the first 5 days for (a) sensitivity experiments on initial conditions, (b) sensitivity experiments to Pr number.

Figure 14 .
Figure 14.The budget of FMSE variance for the first 5 days including: (a) time evolution of the domain mean spatial variance of the vertical-integrated Frozen Moist Static Energy (FMSE) (J 2 /m 4 ) for the sensitivity to Pr number experiments; (b-e)shows the contributions to the FMSE budget by SGS mixing, diabatic feedback (longwave, shortwave and surface flux) and advection (evaluated as a residual term), respectively.Each contribution has an hourly time-step and a 6-hr running average is applied to the horizontal convergence term.

Figure 16 .
Figure 16.The budget of FMSE variance for the first 5 days including: (a) time evolution of the domain mean spatial variance of the vertical-integrated Frozen Moist Static Energy (FMSE) (J 2 /m 4 ) for the horizontal resolution sensitivity experiments;(b-e) shows the contributions to the FMSE budget by SGS mixing, diabatic feedback (longwave, shortwave and surface flux) and advection (evaluated as a residual term), respectively.Each contribution has an hourly time-step and a 6-hr running average is applied to the horizontal convergence term.

Figure 17 .
Figure 17.Sensitivity experiments to horizontal resolution: (a) Vertical profiles of domain-mean cloud fraction avaraged over the first 5 days; (b) Power spectral densities of total water path (TWP) averaged over the whole column and over the first 5 days.The PSD is multiplied by the corresponding wavenumber to have a direct correspondence between the variance of the variables and the areas underneath the curves.The oblique gray dashed line represents the k 5/3 power law.The rightmost vertical gray dashed lines marks the effective resolution of the finer grid resolution model (4 km = 4Δx).The leftmost vertical gray dashed lines marks the effective resolution of the coarser grid resolution model (12 km = 4Δx).

Figure 18 .
Figure 18.Sensitivity experiments to horizontal resolution: (a) average convective distance between updraft over the first 5 days; (b) maximum convective distance between updraft over the first 5 days.Updrafts are defined as grid points where Column Relative Humidity (CRH) exceeds 0.9.Column Relative humidity is defined following the definition by Biagioli and Tompkins (2023) (see their Equation1).A running average with a 24 hr window is used to smooth out the diurnal cycle.

Table 1
Physics Parametrizations Adopted for the WRF and SAM Model Common to All Numerical Experiments as Listed in Table2

Table 2
List of Experiments SILVESTRI ET AL.