Large-Eddy Simulations of Convection Initiation over Heterogeneous, Low Terrain

: Large-eddy simulations are conducted to investigate and physically interpret the impacts of heterogeneous, low terrain on deep-convection initiation (CI). The simulations are based on a case of shallow-to-deep convective transition over the Amazon River basin, and use idealized terrains with varying levels of ruggedness. The terrain is designed by specifying its power-spectral shape in wavenumber space, inverting to physical space assuming random phases for all wave modes, and scaling the terrain to have a peak height of 200 m. For the case in question, these modest terrain ﬁ elds expedite CI by up to 2 – 3 h, largely due to the impacts of the terrain on the size of, and subcloud support for, incipient cumuli. Terrain-induced circulations enhance subcloud kinetic energy on the mesoscale, which is realized as wider and longer-lived subcloud circulations. When the updraft branches of these circulations breach the level of free convection, they initiate wider and more persistent cumuli that subsequently undergo less entrainment-induced cloud dilution and detrainment-induced mass loss. As a result, the clouds become more vigorous and penetrate deeper into the troposphere. Larger-scale terrains are more effective than smaller-scale terrains in promoting CI because they induce larger enhancements in both the width and the persistence of subcloud updrafts.


Introduction
The profound impacts of complex terrain on the initiation and evolution of cumulus convection have received extensive attention, particularly in recent decades.The problem has been addressed using diverse approaches including observational field campaigns (e.g., Raymond and Wilkening 1982;Bougeault et al. 2001;Damiani et al. 2008;Smith et al. 2012), idealized and real-case numerical simulations (e.g., Chu and Lin 2000;Miglietta and Rotunno 2009;Demko and Geerts 2010), and theory (e.g., Crook and Tucker 2005;Kirshbaum and Wang 2014).These efforts have helped to build a reasonable working knowledge of the mechanisms by which orography promotes convection initiation (CI).
As noted in Kirshbaum et al. (2018), orographically induced updrafts capable of initiating deep convection are largely driven by two general mechanisms: (i) impinging flow directly ascending, or detouring around, a given barrier ("mechanical" forcing) and (ii) differential heating over sloping terrain driving thermally direct flow up or down the slopes during the day and night, respectively ("thermal" forcing).These two forcings are not mutually exclusive, and often work in conjunction to enhance low-level ascent (Kirshbaum and Wang 2014).Moreover, terrain features often protrude above nocturnal inversion layers, which locally reduces convective inhibition to facilitate elevated convection.
The amplitude of orographically induced updrafts often, but not always, scales with the terrain relief.In mechanically forced flow, the nondimensional mountain height M 5 Nh m /U is commonly used to assess the flow regime, where N is the Brunt-Väisälä frequency over a relevant layer, U is a representative cross-barrier wind speed, and h m is the mountain relief (Smith 1989).Low-level flow tends to directly ascend the terrain for M , , 1 but split and detour around the terrain for M տ 1.In the former regime, the surface-based updraft amplitude scales with Uh m /a, where a is the terrain width (Kirshbaum and Wang 2014).Analogously, the strength of thermally forced, anabatic flow also tends to scale with h m , due to its control over the terrain-induced baroclinicity (e.g., Crook and Tucker 2005).
Because flow perturbations induced by low, heterogeneous terrain are generally weaker than those over prominent mesoscale ridges (e.g., the Rockies Front Range and the European Alps), the imprint of the latter on climatological flow patterns and precipitation distributions is more dramatic.Thus, the vast majority of studies on orographic convection have focused on larger mountain ranges with reliefs exceeding 1 km.In contrast, relatively few studies have examined the impacts of lower and more random terrain variability on CI.Although such modest terrain variations generate weaker disturbances than their larger-scale counterparts, they may still nontrivially influence the formation and evolution of deep convection, which is highly sensitive to low-level forcing.Moreover, because low-amplitude terrain variations are widespread over land, their collective effects on the global hydrological cycle may be significant.
When smaller-scale terrain heterogeneities have been evaluated in studies of orographic convection, they have largely been considered as perturbations upon larger ridges rather than as interesting features in their own right.For example, multiple studies have analyzed the impacts of small-scale terrain variations on mechanically forced convection over mesoscale mountain ridges.As the impinging flow gradually ascends the mesoscale ridge and approaches saturation, smallscale peaks and ridges can initiate convective cells or bands and anchor them to fixed locations, leading to quasi-stationary convective bands (e.g., Cosma et al. 2002;Kirshbaum et al. 2007a).Similarly, studies on thermally forced convection have argued that valley-ridge systems on mesoscale mountain ranges can concentrate ascent in preferred locations, thus locally focusing the forcing for CI (e.g., Barthlott et al. 2006;Kottmeier et al. 2008).
Terrain is just one of numerous mesoscale features that can promote convection.Another relevant source of topographic variation}heterogeneities in surface cover}has also received extensive attention (e.g., Rozoff et al. 2003;Rieck et al. 2014;Wang et al. 2019).These surface variations, which include land-sea contrasts and spatial variations in land use and/or soil moisture, influence atmospheric boundary layer (ABL) flows in part by modifying the surface energy balance.The resulting mesoscale variations in sensible heating and moistening of the ABL give rise to near-surface buoyancy variations, which, in turn, drive thermal circulations that may effectively initiate deep convection.Spatial variations in surface roughness may also induce low-level convergence zones that force additional ABL ascent.
Although the impacts of land-cover variations on ABL variability and cloud initiation are obvious over flatter land surfaces, they may be obfuscated over complex terrain.Kirshbaum et al. (2016) found that the climatological deep-convection distribution over the southern Mississippi Valley (United States) was mostly controlled by regional terrain variations rather than surface-cover variations.In this region, the maximum terrain relief is approximately 400 m from the center of the valley to the peak of the Ozarks range around 400 km to the west, as the dominant land use varies sharply from mostly cropland within the valley to forest over the mountains.Similarly, the idealized simulations of Imamovic et al. (2017) found that the impacts of regional soil-moisture variations were largely neutralized for 500-m (or taller) mountains.
The objective of this study is to use idealized large-eddy simulations to quantify and interpret the impacts of low, heterogeneous terrain on CI.With the initial conditions, surface fluxes, and mean terrain height held fixed across the experiments, the contributions of different terrain forcing scales on cloud vigor and the timing of CI are systematically examined with the aid of a Lagrangian cloud-core tracking algorithm.Section 2 describes the numerical model and experimental configuration.The model results are presented and physically interpreted in sections 3 and 4, respectively.Section 5 provides the conclusions.

Model setup
Large-eddy simulations are conducted with the cm1 model (Bryan and Fritsch 2002), a fully nonlinear, nonhydrostatic atmospheric model with a terrain-following vertical coordinate and a split time step to stabilize acoustic waves.Time integration is performed with a third-order Runge-Kutta scheme, and advection is fifth-order in the horizontal and vertical, with a fifth-order weighted nonoscillatory (WENO) advection scheme applied to scalars.Subgrid parameterizations include the Morrison two-moment ice microphysics scheme (Morrison et al. 2009) and the 1.5-order prognostic turbulent kinetic energy (TKE) turbulence scheme based on Deardorff (1980).The majority of simulations use a grid spacing of Dx 5Dy 5 250 m and domain size of 120 km 3 120 km in the horizontal, and a stretched vertical grid with a depth of 20 km and a grid spacing that increases from 100 m (over 0-6 km) to 400 m (over 12-20 km) with a smooth transition layer in-between, giving a total of 104 levels.
The simulated flow is based on a case study of tropical convection over the Amazon rain forest, taken from the Large-Scale Biosphere-Atmosphere (LBA) experiment (Rondonia, Brazil) on 23 February 1999.This test case was chosen due to its extensive documentation as well as its prolonged (2-3 h) period of shallow convection prior to CI, the latter of which facilitates analysis of the factors leading to CI. Importantly, however, this case is not chosen to represent the Amazon basin or any other geographic region specifically; it is simply a plausible and interesting test case that serves as a basis for systematic idealized experiments.
In observations and high-resolution simulations of this event, the cloud field deepened slowly, with a growing convective boundary layer (CBL) leading to shallow convection in the morning that transitioned to deep convection at around local noon (e.g., Khairoutdinov and Randall 2006).The model is integrated for 6 h, beginning at 0730 local solar time (LST) and ending at 1330 LST, with model outputs written to file every minute.This high time-resolution is required for the Lagrangian thermal-tracking analysis of section 4b.
The horizontally uniform initial sounding (Fig. 1a), surface fluxes (Fig. 1b), larger-scale radiation tendencies, and momentum relaxation (toward the initial state, with a time scale of 1 h), are adopted from Grabowski et al. (2006).However, unlike Grabowski et al. (2006), who used a free-slip lower surface, we use a "semi-slip" surface where drag is parameterized using a bulk aerodynamic formulation, with a fixed surface roughness length of z 0 5 0.1 m.The inclusion of surface drag improves the realism of the simulated terraininduced circulations, which otherwise may become overamplified (e.g., Richard et al. 1989).Although not shown, sensitivity simulations with a much larger z 0 of 3 m were also performed and gave virtually identical results to the experiments analyzed herein.This insensitivity to z 0 likely stems from the weak low-level winds, the modest terrain heights, and the abovementioned momentum relaxation.The Coriolis force is neglected because the test case is based near the equator.Convective motions are seeded by random initial perturbations in potential temperature on all grid points (maximum amplitude of 0.1 K).The flow is initialized with a horizontally homogeneous passive tracer concentration f that decreases linearly from a value of 0.001 cm 23 at z 5 0 km to zero at the domain top (20 km).
To evaluate the sensitivity of CI to terrain ruggedness, we consider four different terrain fields (h) in the experiments, all of which have a mean height of 100 m: a reference case with a flat surface at 100 m (CTRL) and three complexterrain cases, where h randomly varies from 0 to 200 m over the domain.The complex terrains are designed by specifying their power spectra in Fourier space, assigning random phases to all the modes, inverting back to physical space, and scaling the result to have a maximum relief of 200 m.These spectra all follow slopes of k 1 (where k is the horizontal wavenumber) for scales larger than l c 5 20 km and have zero amplitude at scales below 6D h , where D h is the horizontal grid spacing, to avoid forcing the model at poorly resolved scales.To isolate the impacts of heterogeneous terrain on simulated CI, only the terrain itself is varied between the different simulations.Other aspects of the surface forcing (e.g., surface fluxes, land use, z 0 ) are held fixed across the experiments.
The three complex-terrain cases differ only in the power spectrum of h, or E h (k), between scales of l c and 6D h .Over this wavenumber band, E h (k) follows slopes of k 1 (KP1), to k 25/3 (KM53), to k 23 (KM3), as shown in Fig. 2. The use of E h (k) ∼ k 25/3 as the "central" slope is informed by Kirshbaum et al. (2007b), who found that the Oregon Coastal Range mountains in the western United States followed a nearly k 25/3 spectral slope on the mesoscale, coincidentally similar to the spectral slope of homogeneous, isotropic turbulence.Because the terrains are rescaled to have a maximum relief of 200 m, the three power spectra generally do not overlap in wavenumber space.In physical space, the terrain ruggedness varies greatly between the three spectra, with KP1 dominated by smaller scales (and hence more rugged) and KM3 dominated by larger scales closer to l c (Fig. 3).
Although the magnitude of each Fourier mode is prescribed, the corresponding phase is randomized so that each terrain field represents one of an infinite number of possible realizations for that power spectrum.Given this randomness, one may question whether differences between different cases are owing to genuine physical sensitivities or simply to random variability.To evaluate this question, we run a five-member ensemble of the KM53 case, each identical except for using different random seeds to initialize the Fourier phase spectrum.

Overview of results
The flat-terrain CTRL case undergoes a similar diurnal cloud development to that in Khairoutdinov and Randall (2006), with initiation of shallow cumuli around 2 h, followed by 2-3 h of cloud deepening prior to reaching heights of around 10 km at approximately 5 h (Fig. 4a).Despite the modest terrain relief of only 200 m, the bulk cloud development exhibits notable sensitivities to the terrain.The maximum cloud-top height, after an initial, short-lived jump to around 5 km owing to an ephemeral elevated cloud layer, increases the fastest in the KM3 case, followed by KM53, KP1 and CTRL.At 3 h of model integration, where the differences between the four cases are the largest, the cloud-top height ranges from around 4 km in CTRL to 12 km in KM3.Similarly, the cumuli reach heights of 10 km much earlier (2.8 h) in KM3 than in CTRL (4.5 h).While a similarly strong terrain sensitivity is found for precipitation rate R (Fig. 4c), the sensitivity of volume-integrated cloud mass flux is weaker (Fig. 4b).
Over the 2-3.5-h period where the differences in cloudtop height and R are the largest, the amplitude of the differences between the simulations exceeds the random variability of the KM53 case (as shown by the gray-shaded regions of Fig. 4).While this finding does not offer definitive proof of statistical significance, it nevertheless implies that the variations between the simulations are systematic and not random.Therefore, we conclude that for a given mean terrain height, the presence of terrain complexity generally promotes CI, with the largest effect over larger-scale (KM53 and KM3) terrains.
The above differences in cloud evolution are associated with differences in the low-level flow and cloud water path (CWP), the latter defined as the vertical integral of the cloud liquid and ice water content.At 3 h, the CTRL simulation exhibits a scattered and dense cloud field with generally small CWP (Fig. 5a).While the KP1 simulation also develops a field of scattered cumuli, some cumuli are relatively large with embedded regions of higher CWP than that in CTRL (Fig. 5b).The clouds become more clustered and vigorous in KM53 and KM3, forming longitudinal banded structures with embedded CWP maxima locally exceeding 10 kg m 22 (Figs.5c,d).As will be seen, the largest cumuli tend to initiate on the lee sides of areas of terrain ridges, suggesting the importance of the "leeside convergence" mechanism of Banta (1984).

Physical interpretation
The remainder of this study is devoted to physically interpreting the variations in the timing of CI among the four simulations.Although it is generally expected that complex terrain (however modest) promotes convection initiation, the impacts (if any) of different terrain forcing scales are not immediately obvious.The analysis to follow begins with a brief examination of the sensitivity of boundary layer processes to the terrain heteorgeneities and progresses to an indepth Lagrangian analysis of cloud development.

a. Subcloud processes
To determine the impacts of variations of terrain ruggedness on the bulk convective boundary layer (CBL) properties, we compare profiles of potential temperature (u), water vapor mixing ratio (q y ), zonal (u) and meridional (y) winds, and TKE at 2 h of time integration in Fig. 6.For these analyses, horizontal averages are obtained by first vertically interpolating fields from the native terrain-following coordinates to fixed height levels.When, at a given location, the latter falls below the lowest model grid point, it is assigned a missing value and thereby ignored in the averaging process.Thus, below 200 m, the horizontal averages can be interpreted as conditional averages over atmospheric (rather than subsurface) grid points.Because the surface in the three complex-terrain cases extends above and below the 100-m level, the profiles in these cases extend slightly lower than the uniformly 100-m terrain CTRL case.The TKE is evaluated as the sum of the subgrid TKE (TKE sgs , a prognostic variable in the turbulence closure) and the resolved TKE TKE res 5 u 2 1 y 2 1 w 2 √ =2 , where the primes denote perturbations relative to the corresponding domainwide horizontal average.By 2 h of time integration, the flows have undergone substantial low-level heating and moistening relative to the initial state, along with modest evolution of the low-level winds, below the CBL top at around 400 m (Figs.6a-c).The lowlevel TKE also increases from initially zero to a maximum of around 0.5 m 2 s 22 within the CBL (Fig. 6d).Across the four simulations, the TKE generally increases with the terrain variance, with the smallest values in the flat CTRL case and the largest values in the most rugged KP1 case, with the KM3 and KM53 cases falling in-between.These TKE variations induce commensurate variations in vertical mixing: weaker mixing in CTRL leads to a cooler and moister CBL, while stronger mixing in KP1 leads to a warmer and drier CBL.Although the differences in lower-tropospheric thermodynamic properties between the different simulations are notable (Figs.6a,b), they do not explain the variations in the timing of CI among the simulations.Namely, the two cases with the widest differences in low-level u and q y (CTRL and KP1) are ultimately the most similar in terms of the timing of CI and the resulting precipitation (Fig. 4).Thus, although more rugged terrains exhibit more rapid CBL growth and enhanced mixing between the CBL and free troposphere, differences in low-level vertical flow properties do not directly explain the more rapid CI in KM3 and KM53.Moreover, the variations in low-level winds across the simulations are minimal (Fig. 6c), suggesting that the modest terrain complexity considered herein does not strongly impact the mean lowlevel momentum.
To relate the physical scales of terrain forcing to the subcloud turbulence, we evaluate the power spectra of w LFC , the vertical velocity at the level of free convection (LFC).At a given time, the LFC is determined at a given horizontal grid point (with coordinates x i and y j ) by adiabatically lifting a mean-layer (0-500 m) parcel until it is both saturated and positively buoyant.The local vertical velocity at that height is taken as w LFC (x i , y i ).Discrete 2D Fourier spectra of w LFC (x, y) are taken and multiplied by their complex conjugates, then integrated in rings in horizontal wavenumber k 5 k 2 1 l 2 √ space, where k is the x wavenumber and l is the y wavenumber, to obtain E w (k).These spectra are then averaged over 2-3 h, the period exhibiting the largest intercase differences in cloud-top height (see Fig. 4a).
In all cases, the peak of the E w (k) spectra corresponds to the dominant scale of dry CBL turbulence (l ≈ 1 km), or about twice the mean CBL depth over the 2-3-h period (Fig. 7).The CTRL and KP1 spectra are broadly similar except for the latter having enhanced mesoscale (l . 1 km) energy.Interestingly, this enhancement does not coincide with the peak of the KP1 spectrum at l 5 6Dx (Fig. 2), but rather at larger scales with relatively weak terrain forcing.In the KM53 and KM3 cases, a secondary peak at l 5 20 km emerges due to mesoscale, terrain-driven circulations that develop near the peak terrain forcing scale.At smaller scales (l , 4 km), these spectra are slightly less energetic than CTRL and KP1, indicating slightly weaker CBL turbulence.Thus, while the presence of terrain generally enhances the subcloud forcing in the mesoscale, this enhancement does not cascade down to the microscales.

b. Cloud development
Given the stronger contributions of larger scales to E w (k) in KM53 and KM3 (as just seen in Fig. 7), one may expect the corresponding cloud fields to also be characterized by larger scales.Larger clouds, in turn, tend to be more successful at overcoming the suppressive effects of mixing with dry environmental air as they ascend through the troposphere (e.g., Khairoutdinov and Randall 2006;Rousseau-Rizzi et al. 2017).We thus hypothesize that differences in subcloud forcing, through their impacts on cloud morphology, may hasten the transition from shallow to deep convection.

CLOUD-CORE TRACKING
To evaluate the above hypothesis, we perform a Lagrangian cloud-tracking analysis, where saturated and buoyant thermals (henceforth termed "convective cores," or just "cores") are tracked from their origin near the LFC to their ultimate maximum height.This analysis, described in the appendix, is nearly identical to that described in Rousseau-Rizzi et al. (2017) except that the thresholds to identify core grid points differ slightly, as do the details of the core merging and splitting treatments.To focus on the time period with the largest differences in cloud development (see Fig. 4), we run this analysis only over 2-3 h.
To demonstrate the performance of the tracking scheme, we present snapshots of cloud water path [CWP; defined as the vertical integral of r(q c 1 q i ), where r is the air density] and selected core tracks at three times during the tracking period, over a section of the model domain in the KM3 case (Fig. 8).To minimize clutter, only tracks for cores that reach maximum heights of at least 3 km over 2-3 h are shown.The latter condition eliminates a large number of cores that fail to undergo substantial vertical development.Two of the cores that do reach 3 km, the larger of which grows into the deepest cloud in the entire model domain, initiate in the lee of a small mountain ridge at y ≈ 235 km, after which they widen and intensify as they move southeastward.The individual tracks, which are centered at the location of maximum w within the core, exhibit some chaotic shifts in time, possibly corresponding to the formation of new internal updraft plumes.The largest tracked core also appears to merge with multiple other cores at later times (Fig. 8c), which may help it to widen dramatically over 2-3 h.
We use the core-tracking analysis to evaluate the statistical link between initial core properties and ultimate vertical development.To begin, we present histograms of various core properties, evaluated at the LFC, over the initial phase of development.This phase is defined as the first 10 min of the track, with the analysis restricted to cores based at the local LFC at first detection.For each core, a time series is created over the initial phase or the full core track, whichever is shorter.The initial properties computed from this time series include maximum core area (A LFC ), maximum vertical velocity [(w max ) LFC ], maximum of the mean vertical velocity [(w mean ) LFC ], and time-integrated core-base mass flux (M LFC ).
The total number of cores tracked over 2-3 h is the largest in CTRL (4502) and progressively decreases in KP1 (4090), KM53 (3284), and KM3 (2951).This decrease is primarily associated with a large reduction in the number of smaller cores and a small increase in the number of larger cores (Fig. 9a).The development of larger cores tends to reduce the total cloud population because (i) the volume occupied by a single larger cloud corresponds to that of multiple smaller clouds and (ii) wider clouds have wider descending branches, which inhibit the formation of smaller clouds in between.

FIG. 5.
Comparison of lowest-model-level wind vectors, terrain height (gray fill), and cloud water path (CWP; color fill) for the four simulations, zoomed over a small section of the domain for ease of viewing.
Whereas only modest terrain sensitivities are apparent in the distributions of (w max ) LFC and (w mean ) LFC (Figs. 9b,c), the tails of the distributions of M LFC and A LFC extend to progressively larger values as the terrain transitions from flat (CTRL), to very small-scale and rugged (KP1), to larger-scale (KM53 and KM3) (Figs. 9a,d).Among all of the initial-phase parameters, this trend is the most pronounced for M LFC .Thus, the terrain generally channels more of the CBL ascent into a smaller number of wider cores, particularly in the KM53 and KM3 cases, which increases the averaged intercloud spacing.
The distributions of ultimate maximum heights reached by the cores (H max ), again evaluated over the 2-3-h tracking period (Fig. 10), differ among the four simulations in a similar fashion as the maximum cloud-top height in Fig. 4a.All four distributions peak at H max ≈ 1 km but exhibit a wide variation in maximum core-top heights.While the cores in CTRL fail to reach a height of 4 km, the maximum H max systematically increases in KP1, KM53, and KM3, to around 7.5 km in KM3.These intercase differences in H max are consistent with the corresponding distributions of A LFC and M LFC in Figs.9a and  9d, except for the KM53 and KM3 cases.While the KM53 exhibits the longest tail in A LFC and M LFC , the longest tail of H max is in KM3 case.Possible explanations for this behavior are discussed in section 3c.
Correlations between three initial-phase parameters (viz., A LFC , w max , and M LFC ) and ultimate core-top height (H max ) are presented in Fig. 11.For each initial parameter, five bins are created to cover its simulated range, and distributions of H max within each bin are represented by percentiles (25th, 50th, and 75th).The results show that larger initial values of A LFC , w max , and M LFC are all associated with larger H max , FIG. 6. Low-level, horizontally averaged profiles of (a) u, (b) q y , (c) u and y, and (d) total TKE (resolved plus subgrid) for four simulations at 2 h.
suggesting that these initial parameters may dictate key aspects of the ensuing cloud life cycle.Based on a Spearman correlation analysis (e.g., Wilks 1995), all three of these initialcore properties correlate well with H max , with M LFC exhibiting the strongest correlation (R values ranging from 0.43 to 0.72 over the four simulations), followed by w max (0.28-0.49) and A LFC (0.21-0.42).The corresponding p values are negligible (,0.001), indicating statistically significant correlations at the 99% level for all three initial properties, in all four simulations.Although M LFC exhibits the largest correlation with H max in this analysis, we henceforth use A LFC as our primary initial "control" variable for consistency with past studies (e.g., Dawe and Austin 2012;Rousseau-Rizzi et al. 2017) and its superior predictive skill in the analyses to follow.The positive correlations between initial core properties and ultimate cloud-top height seen above are similar to, but slightly weaker than, those of Rousseau-Rizzi et al. ( 2017) for cumuli initiated by a mesoscale convergence line.These correlations, we hypothesize, reflect a causal mechanism whereby cores with initially larger sizes, vertical velocities, and/or core-base mass fluxes undergo less dilution by mixing with environmental air as they ascend, and thus ultimately become more vigorous and deeper.
This hypothesis is supported by profiles of maximum core buoyancy (b max ; see definition of b in the appendix), vertical velocity (w max ), and total cloud hydrometeor content [(q h ) max ], stratified in bins of A LFC , for the KM3 simulation (Fig. 12).These profiles are calculated by finding the maximum value of the quantity (e.g., b max ) within each tracked core at each vertical level, at each time of the track.The maxima are then conditionally averaged over all cores (and all times) within bins of A LFC .The core-maximum b max , w max , and (q h ) max profiles tend to increase monotonically with A LFC (Fig. 12), with maximum values increasing three-fourfold between the smallest and largest A LFC bins.Similar sensitivities are found if A LFC is replaced by M LFC or w max as the control variable (not shown).
The sensitivities of b max , w max , and (q h ) max to A LFC and M LFC suggests that larger cores with more sustained subcloud updrafts are less diluted by mixing with environmental air than are the smaller cores.This suggestion is evaluated quantitatively by comparing profiles of bulk fractional entrainment and detrainment d, as functions of A LFC .We evaluate and d using the simple "bulk-plume" approach of Betts (1975), using profiles of vertical mass flux (m c ) and a moist-conserved variable f within the cores: averaged over different A LFC bins.Evaluating the resulting profiles for the KM3 case, we find a general decrease in both and d with increasing A LFC (Fig. 13).This trend dominates the profiles below 4 km and continues to hold farther aloft except for a shallow layer where and d are both the lowest in an intermediate A LFC bin (0.6-1.2 3 10 6 m 2 ).This finding suggests that the link between initial-phase parameters and and d is the strongest in the lower to middle troposphere and weakens slightly aloft, the latter of which may stem from the small number of cores that extend above the 4-km level prior to 3 h.

c. Updraft maintenance
The analysis heretofore suggests that heterogeneous terrain promotes CI by generating larger cloud cores, with larger core-base mass fluxes, near the LFC.These larger cores are less suppressed by entrainment and detrainment as they grow, leading to increased midtropospheric cloud vigor and deeper ascent.Although this attractively simple narrative explains the differences between the flat CTRL case and the three terrain cases (KP1, KM53, and KM3), it does not fully explain the sensitivity of CI to the terrain spectrum.For example, whereas the KM53 case has the largest number of large cores during the initial phase (Fig. 9a), CI occurs first in the KM3 case.Moreover, the differences in initial core parameters between the KP1, KM53, and KM3 cases in Fig. 9 are much less prominent than the corresponding differences in H max in Fig. 10.Although one may be tempted to attribute these modest discrepancies to random variability, the narrow spread of the KM53 ensemble in Fig. 4 suggests that the different timing of CI in the three cases is systematic and not random.Thus, other, as-yet unexplained effects must promote faster cloud development over larger-scale terrains.
We hypothesize that, along with initiating larger cloud cores, larger-scale terrains provide a more supportive environment for incipient clouds to develop vertically.First, larger-scale terrains yield a wider spacing between terraininduced updrafts and downdrafts, which prevents newly formed clouds from experiencing terrain-forced descent immediately after forming.Second, the larger-scale terrains generate fewer clouds (and cloud cores) than more rugged terrains, which minimizes competition between neighboring clouds.In more densely packed cloud fields, downdrafts from one cloud may suppress the growth of one or more nearby clouds.Although these two effects are not fully separable in practice, we assume that the former effect manifests as a reduction of core-base mass flux while the latter manifests as enhanced cloud-layer entrainment, evaporation, and detrainment.
Evidence for the first mechanism is provided by time series of mean core-base mass flux (m LFC ) for the KP1, KM53, and KM3 cases (Fig. 14).To focus on the cores with the greatest potential for subsequent growth, only those with A LFC $ 2 3 10 6 m 2 are considered.At each time along the track (t track , which starts at zero for each track), m LFC (t track ) is the averaged LFC mass flux over all tracks satisfying this criterion, excluding those that terminate prior to t track .The resulting time series start at broadly similar values but quickly diverge over time, with KM3 increasing nearly monotonically to large values, KM53 fluctuating with a modest net increase, and KP1 decreasing in time.These results imply that, over larger-scale terrain, larger cores benefit from a more persistent and abundant supply of subcloud air.This may reflect reduced exposure to terrain-induced downdrafts over the core life cycle and/or larger time scales of the subcloud updrafts initiating the cores, both directly related to the increased spatiotemporal forcing scales over larger-scale terrains.
The above argument can be roughly quantified by evaluating the advective time scale of clouds over stationary, terrain-forced updrafts.Consider a simple mechanically forced orographic flow, for which the transit time over such an updraft (t adv ) scales as the ratio of a characteristic forcing scale l c to cloud translation speed U. To estimate U, we consider a shallow cumulus (depth 5 1.5 km) based at the mean LFC over 2-3 h (≈700 m), over which the initial U ≈ 3 m s 21 .The terrain length scale l c is taken as the integral scale of the terrain power spectrum, which is 2.9 km for KP1, 16.5 km for KM53, and 23.4 km for KM3.Using these values, t adv increases from 16 min in KP1 to 130 min in KM3.Despite the great simplicity of this calculation, these time scales are broadly consistent with the experimental trends in Fig. 14, where m LFC in KP1 begins to decrease just 10 min into the track while m LFC in KM3 progressively increases over the 1 h analysis period.
To investigate the second mechanism (cloud-cloud interactions), we apply the bulk-plume entrainment/detrainment analysis to the KP1 simulation (Fig. 15), where the number of tracked cores is 39% larger than in the KM3 case (as shown in Fig. 9).Although the and d profiles again show an inverse sensitivity to A LFC , their values in KP1 are generally larger than those in KM3 (Fig. 13), particularly for the larger A LFC bins.Because these larger cores are the ones that ultimately ascend the deepest (Fig. 11), their increased entrainmentinduced dilution and detrainment-induced mass loss suppresses CI in the KP1 case.Although we defer a detailed analysis of the mechanisms behind these differences in and d to future work, they are consistent with the notion of stronger cloud-cloud interactions, and/or terrain-induced downdrafts extending into the cloud layer, suppressing vertical development in the KP1 case.

d. Robustness tests
Moist convection is a turbulent process that occupies a wide range of spatiotemporal scales.Explicitly resolving this full spectrum of scales is numerically prohibitive, and thus we have adopted a LES approach.This approach is generally acceptable provided that the larger eddies are well resolved and the grid spacing lies well within the inertial subrange, which requires a grid spacing of about 10 times the dominant turbulence scale (e.g., Bryan et al. 2003;Honnert et al. 2011).
In our experiments, this requirement is not met in either the horizontal or vertical directions, where the respective grid spacings (D h 5 250 m and D y 5 100 m) are comparable to the CBL depth (Fig. 6).Thus, it is fair to question whether the simulated CI processes may vary on finer or coarser grids.Moreover, only one initial sounding and surface heating function was considered (Fig. 1), so the conclusions may suffer from a lack of generality.Kirshbaum (2020) systematically studied the sensitivities of both mechanically and thermally forced orographic convection to D h .In flows with well-developed CBLs like the LBA case, D h # 250 m was required for a "converged" solution, wherein further grid refinements did not significantly change the bulk statistics.To evaluate whether this result carries over to the present simulations, we run three tests identical to the FIG.12. Conditionally averaged vertical profiles of (a) b max , (b) w max , and (c) (q h ) max over all cores in various A LFC bins, over 2-3 h in the KM3 simulation.The bin edges are shown in legend.
CTRL, KM53, and KM3 cases except with D h 5 125 m and the terrain refined accordingly.The KP1 case is excluded because, upon refining its terrain, it generated numerical instabilities over very steep slopes with terrain-following coordinates (e.g., Fang and Porté-Agel 2016).Also, to save expense, only one KM53 terrain realization is run at D h 5 125 m, in contrast to the five-member ensemble at D h 5 250 m.
While D h 5 125 m still fails to adequately resolve the larger CBL eddies early in the simulation, it provides a more realistic depiction of turbulent processes in the CBL and above.
Comparing cloud-top heights for these cases to the original simulations with D h 5 250 m (cf.Figs.16a and 4), the trend for CI to occur first in KM3 and last in CTRL is reproduced, with about a 30-min time delay relative to the coarser-gridded original simulations.Similarly, cloud-tracking analysis reveals similar trends as in the original simulations (not shown for brevity).Thus, horizontal grid refinement does not qualitatively change the experimental trends or their interpretations.
We also evaluate the robustness of our findings to ambient winds by rerunning the four cases with the initial wind profile set to zero.Although the winds in the LBA sounding are weak (Fig. 1a), even modest boundary layer winds can have a major impact on topographically forced updrafts (e.g., Kirshbaum and Wang 2014).The resulting bulk time series once again show similar trends in the timing of CI as in the original cases, with the minor exceptions that CI generally occurs slightly earlier and the sensitivity to terrain ruggedness is slightly weakened (Fig. 16b).
The earlier CI in the zero-wind cases is owing to the absence of background vertical shear, which suppresses the vertical development of shallow cumuli by enhancing adverse vertical perturbation pressure gradients (VPPGs) (e.g., Kirshbaum and Straub 2019).The cloud-core-averaged VPPG is approximately halved in the simulations with zero initial wind (not shown), which enhances upward accelerations and facilitates CI.The reduced sensitivity to terrain forcing scale in these cases can be qualitatively explained by the residencetime argument of section 3c.Namely, U 5 0 implies infinite t adv in KP1, KM53, and KM3, which eliminates one factor expediting CI over larger-scale terrains.
Finally, we have run an additional set of three simulations that follow a different approach to varying the terrain ruggedness.Whereas the original simulations varied the slope of the terrain power spectrum, these experiments are based on a single terrain spectrum (KM53, with a k 25/3 slope) and progressively high-pass filter the spectrum to increase the ruggedness.The three cases, KM53-F20, KM53-F10, and KM53-F05, denote filtering of scales larger than 20, 10, and 5 km, respectively.The filtered terrain is then rescaled in physical space to have the same relief (200 m) as the original KM53 case.As the dominant terrain forcing scales become progressively smaller, the cloud development becomes increasingly similar to the KP1 case (Fig. 16c).Thus, these simulations reinforce the notion that smaller-scale/more rugged terrains are less effective than larger-scale terrains at facilitating CI.

Conclusions
Large-eddy simulations of deep-convection initiation (CI) have been performed to determine the impact of heterogeneous, low terrain on cumulus cloud vigor and the timing of CI.The initial conditions and surface heat fluxes were based on an observed case of CI over the Amazon River basin, but the terrain was highly idealized.The latter was designed by specifying a power spectrum in wavenumber space, inverting to physical space using random phases for each wave mode, then scaling the field to have a maximum height of 200 m.Along with a reference flat case (CTRL), three simulations were performed with differing power-spectral slopes in the mesoscale, giving a wide range of topographic ruggedness.
For the case of interest, the terrain notably enhanced cloud vigor and expedited CI by up to 2-3 h, relative to the CTRL case.This promotion of CI was generally owing to the formation of terrain-induced circulations in the subcloud layer that enhanced subcloud kinetic energy on the mesoscale.These broader circulations initiated larger cloud cores (the buoyant and ascending portions of cumuli) at the level of free convection (LFC), leading to less entrainment-induced dilution and detrainment-induced mass loss on their ascent through the troposphere.As a result, the terrain induced the formation of initially wider clouds that become more vigorous and ascended deeper into the troposphere.Although all three terrain fields under consideration promoted CI, larger-scale terrain was more effective at doing so than smaller-scale (i.e., more rugged) terrain.Larger-scale terrain not only initiated wider clouds than did smaller-scale terrain, but it also increased the persistence of the subcloud updrafts feeding the clouds.This increased subcloud support over the cloud life cycle rendered the cloud circulations more plume-like (as opposed to bubble-like).The reduced bulk entrainment rates in the wider, more persistent clouds are consistent with laboratory studies showing smaller dilution rates in updraft plumes than in updraft bubbles (e.g., Blyth 1993).
The above findings were robust to a halving of the horizontal grid spacing D h , a zeroing of the initial winds, and a different method of varying terrain ruggedness, suggesting that they are not specific to the numerical configuration, environmental state, or terrain used herein.However, due to practical limitations, some conclusions may not fully carry over to reality.For one thing, subcloud circulations and microscale turbulence along cloud edges remain under-resolved even on the finest grids considered (D h 5 125 m).While resource limitations precluded further grid refinements, these will be facilitated by future computational advancements.Moreover, the isotropic idealized terrain fields differ from real terrains, which are aperiodic and highly anisotropic.The impacts of this simplification were not quantified but may be significant, and thus will be addressed in future work.Finally, although this study outlined a clear link between terrain scales, CBL updrafts, and associated cloud morphology, it focused primarily on cloud-layer dynamics and did not study the subcloud dynamics in detail.Further analysis of the mechanisms by which heterogeneous terrain modifies the CBL dynamics, across a range of meteorological environments, is thus warranted.Data availability statement.All of the data used for the present analysis are available upon request to the corresponding author.

Lagrangian Tracking of Cloud Cores
The core-tracking algorithm is a postprocessing analysis that relies on overlap in the position of moist thermals between successive 1-min model output times.At a given time, all convective core grid points in the domain are identified by the coincidence of positive vertical velocity w and buoyancy b, along with cloud water plus ice content (q c 1 q i ) exceeding 0.01 g kg 21 .In the preceding, b 5 g u r 2 u r =u r , where the overbar denotes a horizontal average, and q c and q i are the cloud water and cloud ice mixing ratios, respectively.Each interconnected region of core grid points (using 26-point connectivity in 3D) containing 9 or more grid points is treated as an individual core.The 9-point threshold is enforced to avoid very small and ephemeral cores and focus on the more coherent objects.
The above process is repeated 1 min later, and all core objects that overlap with those at the previous minute are treated as continuations of the previous core track, with a unique object identification number.In this fashion, the life cycle of the object is tracked over time.Core mergers occur when, at a given time, a single core overlaps with two or more cores from the previous time.In such cases, the track of the previous core with the most similar volume to the merged core is continued, and the other tracks are terminated.Core splits, in contrast, occur if two or more cores from the current time overlap with a single core from the previous time.Similar to mergers, the split core with the most similar volume to the previous core retains the identification number of the previous core track, while the other cores are treated as new tracks.To ensure that tracked cores originate near cloud base and are not elevated within the cloud layer, we require that core top lie below 3 km for a new core track to be established.

FIG. 1 .
FIG. 1. Configuration of the LBA experiments: (a) skew T-logp plot of the initial sounding and (b) time evolution of prescribed sensible (H) and latent heat (LE) fluxes.In the wind profile in (a), half barbs denote 5 m s 21 and full barbs denote 10 m s 21 .

FIG. 4 .
FIG. 4. Comparison of bulk metrics of convection intensity in the four simulations: (a) maximum cloud-top height, (b) volume-integrated upward cloud mass flux, and (c) domain-averaged precipitation rate.
FIG. 7. Power spectra of w at the LFC averaged over the 2-3-h period, where l is the horizontal wavelength.

FIG. 9 .
FIG. 9. Distributions of properties of tracked convective cores during their initial phase, all evaluated at the local LFC: (a) A LFC , (b) w max , (c) w mean , and (d) M LFC .All quantities are defined in the text.
FIG. 10.Distributions of H max in the four simulations.
FIG. 13.Conditionally averaged vertical profiles of bulk fractional entrainment and detrainment d over all cores in various A LFC bins, over 2-3 h in the KM3 simulation.The bin edges are shown in the legend (units: m 2 ).
Grant NSERC/RGPIN 418372-17 and the U.S. Department of Energy's Atmospheric System Research, an Office of Science Biological and Environmental Research program, under Grant DE-SC0020083.Numerical simulations were performed on the Béluga supercomputer at McGill University, under the auspices of Calcul Québec and Compute Canada.The author is grateful for constructive discussions with Hugh Morrison and John Peters, and thanks three anonymous reviewers for their constructive comments.