On the Sensitivity of Convective Cold Pools to Mesh Resolution

It is well recognized that triggering of convective cells through cold pools (CPs) is key to the organization of convection. Yet, numerous studies have found that both the characterization and parameterization of CP effects in numerical models is cumbersome—in part due to the lack of numerical convergence with respect to the horizontal mesh resolution, Δx, obtained in typical cloud‐resolving simulators. Within a comprehensive numerical convergence study we systematically increase the horizontal resolution in a set of idealized large‐eddy simulations. Our analysis captures key CP processes, namely free propagation, frontal collision and merging of gust fronts. We characterize the numerical convergence of tropospheric moistening rates, gust front vortical strength and propagation speed, and the amplitude of the lobe‐and‐cleft instability. The understanding gained from this analysis may help develop robust subgrid models for CP dynamics.

In an idealized framework-steady, inviscid, and irrotational flow-the propagation speed   can be approximated by   = (2Δ∕) 1∕2 (Benjamin, 1968;von Kármán, 1940), where g is the gravitational acceleration, h is the effective height of the CP, Δρ is the relative density difference between the CP air and that of the surrounding environment, ρ e .However, in order to discuss CP dynamics, the transient behavior between CP initiation and CP dissipation is important, since convective organization depends on the distance traveled by the CP and its ability to trigger new convective rain cells (Falk & van den Heever, 2023).Similarity arguments (Rooney, 2015) predict CP fronts to advance with an approximate square root dependence, where the growth of the mean CP radius  () ∼  1∕2 slows with time.A form where CP radius evolves logarithmically with time was suggested in another recent conceptual model (Romps & Jeevanjee, 2016).
Yet, realistic CP dynamics is more complex: temporally, the CP is affected by the forcing timescale, namely the typical duration of the rainfall event causing the density difference and thus the source of potential energy.Another timescale was recently identified (Meyer & Haerter, 2020) where the nearly laminar-like front would develop strongly turbulent structures.As the GF axisymmetry breaks, the front propagation speed would decrease.Such patterns, which can be characterized as eddies forming along the azimuthal direction along the GF, have previously been referred to as lobe-and-cleft instabilities (Härtel, Carlsson & Thunblom, 2000;Markowski & Richardson, 2010;J. E. Simpson, 1972;Wakimoto, 2001).Based on numerical simulations it was suggested that the activation of these azimuthal features could affect the radial spreading of the CPs by causing a transition from a power law r(t) ∼ t α with α ≈ 0.6 at times before, but α ≈ 0.4 after the transition (Meyer & Haerter, 2020).
Head-on collision effects between CPs have so far not been studied by scaling arguments, likely because the complexity of the fluid dynamics upon collision is prohibitively high.As a consequence, previous studies focused on collisions have resorted to numerical simulations (Feng et al., 2015;Kurowski et al., 2018;Torri & Kuang, 2019).It has recently been found that vertical mass fluxes at the location of collision between two or three CPs are significantly increased and in comprehensive simulations of deep convective cloud fields the locations of collisions were ascribed higher triggering probabilities (Meyer & Haerter, 2020).
Several studies have carefully analyzed the effects of numerical model resolution of CPs (Bryan et al., 2003;Grant & van den Heever, 2016;Hirt et al., 2020;Huang et al., 2018;Moseley et al., 2020;Pressel et al., 2015;Straka et al., 1993;Zuidema et al., 2017), yet the spatial scales involved at the CP gust fronts are often finer than 100 m, as is evident when inspecting the spatial gradient in vertical velocities in observational studies (Kruse et al., 2022).For example, recent results by Grant and van den Heever (2016) advised that simulations should use grid spacings of at most 100 m horizontally and no more than 50 m vertically.Similar conclusions of approximately 100 m spacings were earlier reached by Bryan et al. (2003) and Straka et al. (1993) in their systematic studies of resolutions required to simulate squall lines and density currents, respectively.
In work on horizontally circular CPs spreading in shear-free environments, as we will be discussing them in the current work, azimuthal instabilities can arise in some simulations.It is therefore important to better understand if grid convergence can be reached-thus allowing for a consistent description of energy dissipation within the coherent structures formed at the CP gust fronts.Recent work has suggested that CPs, when more properly resolved, can help inhibit convective self-aggregation (Jeevanjee & Romps, 2013;Muller et al., 2022;Nissen & Haerter, 2021), but may conversely strengthen clustering under diurnal cycle conditions (Haerter et al., 2020;Jensen et al., 2022).
The current study systematically explores resolution effects for CPs generated under idealized initial and boundary conditions, intended to mimic typical realistic rainfall events.The simulation setup is simplified in that it assumes absence of vertical wind-shear-thus allowing for azimuthally symmetric boundary conditions relative 10.1029/2022MS003382 3 of 21 to the origin of the CP.We explore the effects of isolated spreading, head-on collisions as well as CPs merging, thus addressing processes and effects relevant to cloud organization under a range of typical atmospheric conditions.
Our key findings are: 1. So-called "cloud resolving" grid resolutions,  (1 ) are far from sufficient in the description of the CP gust front including its vortical dynamics; 2. The azimuthal mode is activated for Δx values near and finer than  (100 ) and azimuthal and turbulent processes increasingly cause dissipation and thus reduction in kinetic energy at finer scales; 3. Merging of CPs, for example, within MCS, can bypass the slowdown of radial spreading, and a fast spread with r(t) ∼ t is nearly reached for merging CPs-thus allowing MCS to spread further into their environmentand the flow in their wake subject to a better-sustained mechanical lifting caused by trapped reflection waves.

Large-Eddy Simulation
This study builds on the recent work by Meyer and Haerter (2020) and uses the same numerical methods, summarized below.The simulations use the large-eddy simulation (LES) model PyCLES developed by Pressel et al. (2015).The model solves an anelastic form of the momentum and moist entropy conservative equations on an orthogonal mesh.A 5th-order weighted essentially non-oscillatory (WENO) scheme (Jiang & Shu, 1996) is used, which keeps grid-scale oscillations low even in the vicinity of discontinuities, such as the front of the CPs, at the cost of higher dissipation (Pressel et al., 2017).The time stepping is determined dynamically using a fixed Courant-Friedrichs-Lewy (CFL) number, ensuring stability of the explicit time stepping scheme.The sub-grid scale turbulent closure uses a Smagorinsky scheme (Smagorinsky, 1963).Surface fluxes of momentum, specific humidity and entropy are modeled by exchange laws as in Pressel et al. (2015) and evaluated at the lowest model layer using Monin-Obukhov similarity theory (Byun, 1990).Exchange coefficients are chosen in accordance with simulations of warm rain convective clouds (c h = 0.001094 for heat, c q = 0.001133 for moisture, and c m = 0.001229 for momentum flux;vanZanten et al., 2011).Importantly, this means a no-slip surface boundary condition is only implicitly applied to the flow by means of an increased inner-cell viscous dissipative flux.The top region between 3 and 4 km altitude is treated as a sponge layer meant to prevent any spurious reflections from disturbing the inner solution.The data is exported every Δt = 30 s, which we herein conveniently refer to as the timestep.

Setup of the Simulations
In order to study the effects of mesh resolution on single CP spread, collision and merging, we consider the following three types of simulations: 1. Free propagation: A 24 × 24 km 2 horizontally periodic domain is constructed with a cylindrical forcing region placed at its center.The forcing is prescribed during a time interval of τ = 900 s, following the procedure described in Section 2.3.This setup results in a column of cold and moist air at the domain center which, upon collapsing onto the surface, creates a CP that spreads freely within an approximate 10 km radius over the 1-hr runtime.

Frontal collision:
To simulate the frontal collision between two CP fronts, we adapt the horizontally periodic domain to 12 × 20 km 2 and center the cylindrical forcing region at (x, y) = (0, 10) km.Owing to the domain's periodicity, this splits the cooling region into two half-cylinders.The domain symmetry causes the collapsing air to meet half-way through the domain.These simulations are also run for 1 hr in physical time.3. Merging: A narrow 96 × 16 km 2 horizontally periodic domain with a (τ = 3,600 s)-long cylindrical forcing placed at its center.This longer forcing time provides a greater potential energy to the CP, which sustains its propagation further away than for the two other cases and allows us to capture the incidental merging between fronts.These simulations are run for 4 hr in physical time.Additionally, a reference single-cold pool 96 × 96 km 2 case is run at 100 m-resolution to compare with the merged front configuration.
All simulations are initialized in an idle state, with an atmospheric stratification profile as adopted from Grant and van den Heever ( 2018 (1) Importantly, the atmosphere is sufficiently dry such that no condensation will occur in any of our simulations and q t = q v , where q t is the total specific humidity, and Finally, all simulations use a structured orthogonal mesh with the same vertical resolution Δz = 25 m and explore a range of horizontal resolutions Δx ∈ [25,50,100,200,400,800] m.A domain height L z of 4 km is chosen after performing a sensitivity study demonstrating this domain height to be sufficient to prevent excessive strain on the lower-level circulation.It is reminded that the domain is topped by a sponge layer, from 3 to 4 km altitude.
Overall, this setup allows us to study the three aforementioned processes in a numerical environment, shielded from external perturbations or large-scale motion, in an effort to expand our understanding of CP dynamics.This understanding is, however, constrained by the unique atmospheric conditions we chose to study.This practical choice was motivated by the study's computational cost.All runs are summarized in Table 1.

Cold Pool Forcing
The forcing consists of a cooling and moistening that is applied as fixed tendency terms over a cylindrical volume within the computational domain, labeled V CP , and a certain period of time τ.V CP extends from the surface to the top of the neutrally stratified planetary boundary layer (z* = 1 km) at a radius of r* = 2 km.
The cooling is applied as a diabatic tendency term to the balance equation of entropy where the dot denotes a time-derivative, δT is the cooling rate in kelvins per second and c p,m (q t ) = (1 − q t )c p,d + q t c p,v is the specific heat of moist air at constant pressure.
To mimic how the CP is generated by the evaporation of rain, we evaluate how much liquid water has to evaporate to produce the cooling rate δT.To this end, an equivalent moistening term δq t is computed.Note that the cooling is applied homogeneously through V CP , while the moistening rate varies with height due to the dependence of the latent heat L on the temperature T, which in turn decreases with height.For the sake of simplicity, we neglect drag by rain drops throughout this study.
The equivalent moistening rate is evaluated as following.For a given volume V of air, the water vapor specific humidity is defined as q v = m v /m t = ρ v /ρ t , with m v and m t the specific mass of water vapor and total air, respectively, and ρ v and ρ t the corresponding mass densities.A small change in moisture, δq v , can then be expressed as: with m d the mass of dry air, such that m t = m d + m v .Consequentially, it is: For an air parcel of volume V, we can now relate the amount of cooling δT to the corresponding mass of evaporated water δm v (5) with the latent heat of vapourization, L, a function of temperature T, and the heat capacity of moist air c p,m (q t ).
Hence, the change in the water vapor specific humidity is where the latter equality holds for sub-saturated air, as is the case throughout this study (condensation never occurs by design).
On the level of prognostic variables, this translates into an imposed tendency of the specific entropy s and the total specific humidity q t , namely: Here, s v is the specific entropy of water vapor, a function of the partial pressure of water vapor p v , and s d the specific entropy of dry air, a function of the partial pressure of dry air p d .
In this study, we arbitrarily set δT to −17 K hr −1 to replicate a similar temperature anomaly after the cooling period of τ = 15 min as in Meyer and Haerter (2020).Corrections are applied to compensate for the changing volume of the cooling cylinder volume at varying grid resolution, such that the added potential energy remains rigorously constant over all resolutions.Considering the resulting maximum temperature depression of approx.−4.25 K, the simulated CPs are colder (and hence stronger) than oceanic CPs (Addis et al., 1984;Zuidema et al., 2012Zuidema et al., , 2017)), but correspond to typical CPs from moderate single convective cells over the (mid-latitudinal) continent (Addis et al., 2021;Kruse et al., 2022).

Cold Pool Formation and Propagation
The first process to be scrutinized in this study is that of CP propagation, as described in Section 2.2.A qualitative view of the simulation is provided in Figure 1 at several key timesteps.An animation is provided as Movie S1.The CP is visualized from the potential temperature anomaly using the θ′ = −0.1 K isosurfaces, colored by altitude.This threshold of −0.1 K was arbitrarily chosen to capture the CP edges while filtering out the minor random perturbations used during initialization.Conversely, a higher threshold would only reveal the interior of the CP and not provide more insight.At t = 300 s, the cylindrical forcing zone is observed at the center of the domain, forming a cold column of air reaching 1 km altitude.At t = 900 s, the base of the column has started to collapse, yielding a circular front spreading outwards.Note that the forcing is still operating: the cylindrical cooling zone continues to reach the 1 km-level.At t = 1,500 s, the forcing has stopped and the central column of cold and moist air has rapidly collapsed.The front spreads further away from the epicenter and the GF exhibits azimuthal perturbations in both vertical and radial directions.These correspond to lobe-and-cleft (LC) instabilities, which are the dominant 3D instabilities occurring in gravity currents (Cantero et al., 2005;J. E. Simpson, 1982).At t = 2,100 s, the CP continues spreading outward, with a noticeable growth of the LC modes.Note the regularly-spaced red-colored protuberances indicating large vertical deformations of the leading front.By t = 2,700 s, the central column of moist air has totally vacated its initial zone, and been replaced by entrained dry air.This has resulted in a central θ′ > −0.1 K region (in fact, a warm anomaly) as indicated by the isosurface's disappearance.Finally, at t = 3,300 s, the front reaches the edges of the computational domain.By then, it has already started self-interacting owing to the domain's periodicity.This marks the end of the free-propagation phase.As a result, the remainder of the study only considers simulation time up to t = 3,000 s when the front is still away from the border.Importantly, we detected no other kind of spurious reflections at the domain's edge during the timeframe of interest, such as density-waves.

Cold Pool Visualization
A first qualitative view of the CP's sensitivity to mesh resolution is presented in Figure 2. The near-surface (z = 12.5 m) vertical velocity w contours are plotted at various timesteps.As the mesh coarsens, numerical dissipation increases, as is expected from a dissipative scheme such as 5th-order WENO (Balsara & Shu, 2000).This broadens the front and reduces the velocity amplitudes.Additionally, the increase in horizontal mesh element size appears to cancel out the growth of azimuthal asymmetries: the Δx ∈ [200, 400, 800] (lower half) contours retain azimuthal symmetry (or flow laminarity) until the front reaches the domain's boundaries.This proves that the activation of the LC instability is physically-driven by wavelengths of size strictly lower than 400 m (i.e., the 200 m-mesh Nyquist cut-off).Further, the high-resolution simulations present ring-like structures trailing in the  wake of the leading front.These correspond to secondary vortices detaching from the main leading vortex as a result of the continuous forcing until τ = 900 s.They do not form in the Δx > 100 m cases, where the flow in the front's wake appears quiescent.Hence, while numerical dissipation erodes the fronts in the Δx > 100 m cases, the absence of secondary vortex and LC instabilities helps those fronts remain more coherent, relative to those in the Δx ≤ 100 m cases.

Development of the Lobe-and-Cleft Instability
To analyze the LC instability, the CP frontline is extracted from the surface level and plotted for every single timestep in Figure 3, still for the Δx = 25 m resolution.The resulting flower-like structure reveals the front's progressive deformation over the CP lifetime, which can be divided into two stages.At first, the CP spreads in a remarkably axi-symmetrical fashion.That is, its local propagation speed is invariant along its azimuthal angle ϕ which results in a series of perfectly concentric circles centered around the forcing region.Importantly, this stage corresponds to the accelerating stage as is indicated both by the red-shifting (see color scheme) and increasing distance between neighboring circles.The black circle marks the t = 900 s time step, when the forcing is turned off.This coincidentally marks the end of the accelerating, laminar-like phase and the transition onto the second stage.This second stage starts with the onset of the LC instability which grows rapidly, causing a deformation of the frontline and its rapid deceleration.The lobes correspond to regions of high velocities with respect to the mean front propagation speed: they are ahead of the mean front radial location,  ( ) > () , where r(ϕ, t) is defined as the front's radial distance from the domain's center and   denotes the azimuthal average.Conversely, the clefts are the regions suffering from the strongest deceleration, visualized as dark blue wrinkles.This imbalance creates pairs of counter-rotating streamwise vortices which have positive (respectively negative) vertical mass fluxes around the lobes (respectively clefts), as was visualized by Dai and Huang (2022).We observe that the clefts can merge with one another but never disappear (see enlarged sub-panels in Figure 3).We interpret this as a direct consequence of the apparently irreversible local flow slow-down that caused them in the first place.In absence of a forcing that would re-accelerate the flow, the imbalance cannot disappear.Conversely, the lobes grow over time as their vortex pair lifts them above the front where they can freely expand.
Although the LC instability has long been studied in the context of density current (Cantero et al., 2005;Dai & Huang, 2022;Härtel, Carlsson & Thunblom, 2000;Härtel, Meiburg & Necker, 2000 et al., 2022), there exists no consensus as to its triggering mechanism.The most-recognized explanation for the formation of LC waves attributes it to the no-slip surface condition: the shear stress causes the bulk of the gravity current to be raised some distance above the slower near-surface front (J.E. Simpson, 1972Simpson, , 1982)).The cold front head is therefore lifted over a small parcel of the surrounding air which triggers a buoyancy-driven LC instability.This theory has gained much traction thanks to mechanism-denial experiments carried out by means of direct numerical simulations (Cantero et al., 2005;Härtel, Carlsson & Thunblom, 2000), showing that a slip-wall condition disabled the LC instability.More recent numerical studies further support this theory (Dai & Huang, 2022;Xie et al., 2019), although it must be reminded that these expensive numerical studies were limited to low Reynolds numbers differing from real-world CPs by many orders of magnitude.Conversely, Horner-Devine and Chickadel (2017) recently observed LC instabilities forming in absence of strong shear stress, giving weight to an alternative theory from Parsons (1998): the LC mode can form as a response of the breakdown of a Kelvin-Helmholtz (KH) wave in the upper layer of the gravity current.As the KH billows collapse onto the front, they initiate streamwise vortices growing into the LC instability.This mechanism would be dominant at high Reynolds numbers (Horner-Devine & Chickadel, 2017).Interestingly, KH waves appear to precede the onset of the LC mode in Xie et al. (2019) as well, which complicates the determination of its forming mechanism.Determining the origin of the LC instability is beyond the scope of this study, but we make the following observations.First, while the Monin-Obhukov model artificially increases dissipation to represent surface drag, an underlying slip surface boundary conditions is used in essence.Indeed, we observe no prominent front over-head.Yet, an LC instability unquestionably develops at high resolution.Second, the Reynolds number Re f based on the front thickness and velocity at simulation half-time t = 1,500 s, after the LC mode has activated, is evaluated as Re f ≈ 3 • 10 8 .In light of these observations, the existence of a high-Reynolds secondary mechanism (Parsons, 1998) appears plausible.
As a result of the LC instability and growth of streamwise vorticity, momentum is transferred from the radial into the azimuthal direction.This is illustrated in Figure 4a which presents the frontlines projected in an azimuthal coordinate system-(Δr, Δϕ) = (20 m, 1°)-and colored by the local azimuthal velocities u ϕ evaluated by linear interpolation.Clefts are easily identified in the form of frontline wrinkles, that is, low-velocity regions trailing behind the frontline.Note the change of sign from positive to negative around the clefts indicating a convergence area.Conversely, the lobes constitute the prominent u ϕ = 0 m/s regions marked by a negative-to-positive u ϕ shift indicating a diverging flow.As was observed in Figure 2, the growth rate of the LC mode varies enormously between the four most refined cases presented in Figure 4a.For instance, the first clefts in the Δx = 25 m case appear upstream of those for the Δx = 100 m case.
In an effort to better quantify the LC length scale, we resort to Fourier analysis.A 2π-periodic azimuthal signal for u ϕ is extracted from the last t = 3,000 s frontline.Its single-sided amplitude spectrum is plotted in Figure 4b for all spatial resolutions.As was observed in Figure 2, the azimuthal mode is much less energized as the mesh coarsens.More surprisingly, however, the finer-resolution amplitude spectra appear quite broadband, while we expected strong peaks from the onset of the LC instability.This suggests that the lobes exist at a variety of scales.We mark with double-arrows what we interpret as the most-energized bands, and observe that there is no clear convergence as spatial resolution improves.The mid-to-high wavelengths seem, however, well-resolved by Δx = 100 m.
The time of activation can also be evaluated from the frontline signals (Figure 4c).To this end, we use the standard deviation of the CP radius r along the azimuthal direction, σ r(t,ϕ) , as measure of the symmetry-breaking caused by the LC instability.By normalizing this value to the mesh resolution, it is possible to accurately define the time of activation which increases with mesh element size.It is reminded that the forcing ends at t = 900 s, which seemingly coincides with the onset of the LC mode.To test whether the end of the forcing impacts the LC mode, a 1h-long forcing case is run at 100 m resolution (see dashed line in Figure 4c).It appears that the onset of the LC instability occurs roughly at the same time as for the shorter forcing period, although its growth is impeded by the continuous radially-oriented forcing.The transition time τ LC , denoted as when the indicator cross an arbitrary threshold of 1.5, is marked for the finest cases.It increases linearly as resolution improves, which suggests that convergence has not been reached.It seems to approximately evolve as a linear function of Δx, which suggests that the 200 m-resolution case might have ended up developing an LC instability had the domain been larger.
Finally, the frontline propagation speed, U r , is evaluated at every timestep as the ϕ-averaged radius rate-of-change.It is plotted for all cases in Figure 4d.The initial discrepancies between cases is an artifact of the mesh coarsening, as the frontline is poorly approximated by squared elements.This spurious effect becomes less pronounced as the CP spreads.It then becomes clear that there is very little difference in advective group velocities, indicating that the GF speed is resilient to mesh coarsening.Consequently, the LC mode has little-to-no impact on the front propagation speed.This contradicts the findings of Meyer and Haerter (2020), although their configuration lacked a constant forcing and initialized a column of cold air to trigger a density-current.The change of front velocity could have been coincidental with the appearance of the LC instability.Importantly, three distinct phases are identified from Figure 4d, just as described in Yuan et al. ( 2022): 1.An acceleration over the first 250 s, akin to the slumping phase from Yuan et al. (2022).2. A constant-acceleration plateau of a few hundred seconds, starting shortly after the end of the forcing.It would correspond to the inertia phase from Yuan et al. (2022).3. A deceleration, which should be followed by a collapse and of the GF, leading to its eventual dissipation.This would occur outside the bound of the computational domain and is not captured here.It corresponds to the viscous phase as defined in Yuan et al. (2022).

Evolution of Gust Front Vorticity
A key feature of the CP gust front is its ability to lift the air it passes through.This lifting mechanism is associated with the strong horizontal vortical activity embedded within the front.Hence, the evaluation of the leading vortex's size and strength is a relevant metric to assess the ability of a mesh to satisfactorily resolve a CP.
The azimuthally averaged r-z composite in Figure 5 shows the radially advancing gust front (GF).The GF is composed of cold and dense air, which spreads along the surface under gravity.All simulations of different mesh resolutions show a vortex embedded within the GF, by which air is lifted mechanically upon its passage.This primary vortex induces a positive water mass anomaly above 1 km under the assumption that the boundary layer is moister than the highest levels.Interestingly, the streamlines show how the vertical extent of the vortex ring exceeds 2 km.As is expected from vorticity conservation, a counter-rotating vortex trails behind the front, which in turn reinforces the downward flow within the CP wake which closes the contour of positive moisture anomaly.Importantly, this secondary vortex is not associated with a strong positive water mass perturbation as it first mechanically lowers the incident air before lifting it up.The main effect of decreased resolution is to deteriorate the internal structure of the primary vortex ring, with positive water mass anomalies at the highest resolution markedly diminishing and smearing out horizontally as the coarsest resolution is approached.Note the red "+" symbols, marking the vortex ring's maximum and minimal vertical velocity locations, moving closer to each other between Δx = 25 m and Δx = 100 m and subsequently further away from another.
It is important to numerically resolve the CP gust front width by a sufficient number of discretized grid boxes, with literature (Marburg, 2002) recommending at least six elements per wavelength for resolving one-dimensional wave transport with a high-order discretization scheme such as the one used here.Indeed, at coarse spatial resolution, Δx ≥ 200 m, this requirement it not met (Figure 6a), where the GF is only covered by between 2.5 and 4.3 elements (meshpoints).The increasing resolution between Δx = 800 m and Δx = 200 m physically manifests itself in a spatially more confined, yet moister, vortex ring (Figure 5).This lifts more moist air into the free troposphere, as quantified in Figure 6b by plotting G, defined as the maximum column-integrated moisture gain-caused by this bubble-averaged over the simulation time, G.It is calculated as: where   is the time-window of interest starting from when the forcing ends (t = 900 s) and ending as the front nears the domain edge (t = 3,000 s).
At Δx ≈ 100 m the LC mode activates, and the associated symmetry breaking along the frontline gives rise to a range of secondary processes, such as trailing vortices that even detach from the primary vortex (compare: Figure 5, Δx = 25 m, black arrows).We interpret the effect of the LC instability, which effectively transfers energy to an azimuthal mode, as weakening the kinetic energy of the circulation in the r-z plane.This would explain the increased moisture gain G for the intermediary resolutions.This interpretation is further supported by considering the integrated vortical strength (Figure 6c), a quantity that decreases with mesh resolution after the onset of the LC instability (t ≥ 1,800 s, Δx ≤ 100 m).This quantity corresponds to the spatially-integrated rigid-body rotational rate ω r , a quantity extracted from the velocity field, akin to vorticity.However, contrary to vorticity, it does not account for shear and strain rates: one can therefore quantify the pure rotational strength in a stratified, shear-dominated flow such as this one.The reader is referred to Kolář (2007) for a presentation of the triple decomposition method.In summary, in the later stage of expansion (t ≥ 1,800 s), CPs simulated at too low resolution may either overestimate (intermediate resolution Δx ≈ 100 m) or underestimate (very low resolution) vorticity relative to the Δx = 25 m case.During the early stage (t < 1,800 s) before the LC mode activates, vorticity increases monotically with mesh resolution as numerical dissipation decreases.

Scalability of the Upward Mass Water Fluxes
CPs constitute agents that can mediate interactions between co-occurring deep convective rain cells-thus allowing existing convective cells to form new ones in their surroundings.Since new cells are often triggered by moisture lifted vertically, it is important to investigate the upward moisture fluxes brought about by the CPs as they spread (Figure 7).The comparison of the probability density functions (PDFs) of upward moisture mass flux for different resolutions shows that they essentially coincide for low values of upward moisture mass flux, whereas they depart for higher moisture fluxes.Using the highest-resolution PDF (Δx = 25 m) as a reference, it can be seen that these departures are systematic and predictable as one moves from resolu tion to resolution.Before the points of departure, all distribution functions have an estimated power law exponent of ≈−2.4.For an unbounded power law distribution of this exponent extreme fluctuations in the higher moments of the distribution, including the second moment and thus the variance, 〈J 2 〉 − 〈J〉 2 , would be implied.With a divergent variance, extreme values of J would be  statistically unpredictable.All distributions for the different spatial resolutions Δx do show bounds, that is, a scale break of the power law distribution, yet, as Δx → 0, this bound seems to be systematically pushed toward larger and larger J. Quantitatively, the departure of the lowest resolution curve (Δx = 800 m) occurs at moisture fluxes that are more than an order of magnitude smaller than those still resolved at Δx = 25 m.This implies that especially extreme moisture fluctuations, brought about by CP spreading, may systematically be cut off in low-resolution simulations.Notably, the decay of cut-off fluxes with horizontal resolution is remarkably linear: compared to the highest simulated resolution Δx = 25 m, the Δx = 50 m fluxes are resolved up to the 99.6th percentile (only the top 0.4% was dissipated) whereas the Δx = 800 m fluxes are resolved up to the 81.0th percentile.The other resolution fall along the line formed by these two limiting cases.Considering that new rain cells may in particular be triggered by the larger moisture fluxes, this reduction may have substantial consequences for realistic simulation of the cloud-to-cloud interaction dynamics.Finally, these results suggest that (a) numerical convergence of the bulk vertical mass flow rate is attained for Δx ≤ 100 m (99th-percentile resolved), (b) the tail of the distributions, corresponding to the unresolved fluxes, follows a power law which could serve as the basis for a CP parameterization, akin to turbulence closure scheme.

Frontal Collision Between Gust Fronts
CPs spreading freely in space are rarely found in realistic simulations of the convective cloud field (Haerter et al., 2019;Nissen & Haerter, 2021).Rather, nearly all CPs collide with other CPs at some stage during their expansion.Thus, we now simulate idealized collisions between two identical CPs, which only differ as a result of the small random perturbations imposed during initialization for symmetry-breaking purposes (Figure 8).
In our simulation setup the horizontally cyclic boundary conditions cause a given CP to interact with itself, since a collision of the CP's GF occurs across the boundary.Based on the approximate mean distance between adjacent cells in the simulations of Haerter et al. (2019), we use a domain width of 12 km.The qualitative finding is that, at locations of collisions, the perturbation in potential temperature extends to a greater height than at locations where the GF spreads freely.Geometrically, the initial pattern formed by collisions is that of a straight line, along which perturbations are similar (near t = 1,500 to 1,800 s), akin to those within Voronoi diagrams (Haerter et al., 2019;Kim et al., 2001).However, at later stages during the collision (e.g., t = 2,100 s), the strongest perturbations are located within a point-like geometry, that is, the two points defined by the intersecting circular gust fronts.Our interpretation of this finding is that the points located at the "edge" of collision constitute a one-dimensional "gust front," propagating along the line formed by the collision.For colliding circles of equal radius r, which spread radially at velocity v r , with centers separated by a distance 2d, the velocity at which the collision line grows is given by , which can be checked by noting that  () 2 =  2 +  2 .Thus, these experience larger spreading velocities than do other portions of the GF.In the later stages of the collision the density current is reflected back and gives the visual impression of the initial circles simply expanding further.Yet, seeding tracer particles confirms that there is essentially no mass transport across the line of collision (not shown), and there is therefore no mixing between the two CPs.
In order to assess the capability of collisions of enabling new convective cells, we characterize the upward moisture perturbations along two pathways (Figure 9).The so-called collision path corresponds to the y = 10-km plane, and is centered on the shortest time-to-collision pathway.Conversely, the free path corresponds to the x = 0-km plane, and is centered on the longest time-to-collision pathway.By the time the collision has been resolved on the former path, the later path would still be collision-free.
First, we observe that the peak in free-tropospheric water anomaly occurs prior to the collision of the two vortex rings, as seen in the   ′  contour levels (Figure 9a, t = 1,200 s).Thus, new convective cells might sometimes be triggered ahead of gust fronts collisions, rather than after them.This is further quantified by evaluating the quantity of free-tropospheric water mass perturbation, Q v , in Figure 9b.Q v is calculated as the integral of water mass perturbation in a volume zone of size (ΔX, ΔY, ΔZ) = (3,200, 1,600, z > 1,000) m that is placed at 5 km from the epicenter and is perpendicular to both pathways.This 5 km corresponds to the exact center of the collision path, hence is centered on the collision, while it only captures one front propagating along the free path.Q v is found to grow about twice faster along the collision path compared to the freely propagating CP, due to both fronts simultaneously entering the zone of interest around t = 400 s.Overall, the maximum perturbation is about twice as large in the collision case, which indicates that collisions are not 100% efficient in converting horizontal momentum into vertical displacements.It is encouraging that this crucial metric is similarly evaluated by all but the lowest resolution, where a notable reduction is visible.

Incidental Merging Between Cold Pools-Gust Front Speed
Mesoscale convective systems may form when multiple deep convective cells merge to generate a joint, quickly-moving GF, to be built up (Haerter et al., 2020;James et al., 2006;Skamarock et al., 1994).A more recent followup work adds to this in suggesting that this effect becomes more pronounced as horizontal resolution increases (Jensen et al., 2022).As an additional effect, the formation of a macro-CP by merging of several CPs was seen to cause persistent, multi-day, drying over large regions (Jensen et al., 2022).The resultant combined CP covers substantially larger horizontal areas, reaches greater height, and is often longer lived than individual CPs (Feng et al., 2015).Crucially, due to its greater height and spreading speed, the combined CP may lift moist air masses to the level of free convection without requiring collisions.Using the idealized elongated geometry described in Section 2.2, we address how combined CPs form and spread.The merging processes is presented in Figure 10.
The initial state is similar to the one used before (Figure 10a compared to Figure 8).The raincell epicenter is located at the origin coordinates,   ≡ (0, 0) .We refer to the first point of collision as A e (0) ≡ (0, d)-subscript e for edge- where = 8 km defines the distance between the epicenter and the domain's edge, such that   = |(0) − | .After the collision occurs, the front at A e (t) quickly starts "catching up" with the central position, A c -subscript c for central-which spreads freely in a circular fashion.Whereas this "catching up" is expected on geometrical grounds, it is interesting to note that the expansion in y-direction along the domain edge actually overtakes that of the domain symmetry axis.We attribute this to the symmetrical periodic boundary conditions which converts all the momentum in the y-direction into the vertical z and axial x directions, thereby accelerating the x-velocity and increasing density anomalies near the edges, which further reinforces the density current.As a result, A e grows faster than A c , and overcomes it.Eventually, near t = 10,800 s, A c (t) is similarly reinforced by the central collision of the reflected waves, and its x-coordinate manages to surpass that of A e .Thus, an oscillatory internal front dynamics evolves.
Again turning to a quantitative analysis, we now measure the front's propagation in terms of its mean y-coordinate and plot this position, termed δ(t), for the various mesh resolutions versus time (Figure 11a).
Before the collision, the merged front simulations all follow the same front dynamics as the freely-propagating simulation (Figure 11a).This is unsurprising as the two configurations are, up to this point, identical, bearing the negligible impact of virtually-indistinguishable density waves, owing to the domain's periodicity.
After the end of the forcing at t = 3,600 s, the single and merged front decelerate.In contrast to the free propagation case (Figure 1), the impact of resolution on δ(t) eventually becomes noticeable: propagation systematically slows as resolution coarsens with final values of δ(4h) ranging from 36.4 to 39.3 km for Δx = 800 and 25 m,  respectively.Yet, even the coarsest merged front reaches further than the reference case (Δx = 100 m) for the CP configuration after the whole 4h-long simulation (Figure 10a, inset).
Finally, we compare the single and merged fronts propagation rates-both with Δx = 100-once the forcing has ended, in Figure 4d.The front locations appear to be initially well resolved by power laws of exponents γ = 0.74 and β = 0.83, for the single and merged fronts, respectively.These exponents are fitted to the solution at δ(t = 3,600 s) (dashed lines).In both cases, the front's deceleration rate appears to increase toward the end of the simulation, when the solid and dashed lines separate around t = 7,800 and 11,000 s for the single and merged fronts, respectively.We deduce from the higher exponent and longer holding time that the merging of neighboring CP, caused by the channel-like configuration, has reinforced the central front, allowing it to propagate further away and faster.

Incidental Merging Between Cold Pools-Mechanical Lifting
After assessing the merged front traveling speed as a function of mesh resolution, we now turn our attention to the merged front's strength.That is, we study the merging's capability to enhance upward moisture fluxes, just as was done with the single CP case (Figure 7).This analysis is presented in Figure 12a in similar fashion.Strikingly, the PDF for all resolutions but the coarsest one share a self-similar profile characterized by the same −2.4power-law.This is particularly noteworthy as it implies any flux-enhanced model based on this power law would be equally suited to resolve isolated and merged CPs alike.As both processes are bound to occur in a realistic simulation of the atmosphere with multiple CPs typically erupting quasi-simultaneously, this is reassuring.Another point of interest to modelers is that the cut-off flux appears, once again, predictable as function of resolution.This time however, the power-law scaling is of lower magnitude: hence, merged fronts suffer less from grid coarsening than do isolated CPs.
Finally, we explore the amplitude of the mechanical lifting for all resolution and the 96 × 96 km 2 single-CP reference case.To do so, we define a scalar metric in space and time corresponding to the 90th percentile of front altitude z 90 .This percentile is evaluated amongst the z-population in all radially-extending bins of size 800 m starting from the epicenter.We define an equivalent radial coordinate r* for the channel configuration which simply consists of a cylindrical coordinate transform around the Cartesian origin.For the 96 × 96-km 2 configuration, r = r*.The fronts are those shown in Figure 10, that is, defined as the θ′ = −0.1 K isosurfaces.
The profiles are shown in Figure 13 with the dashed line corresponding to the reference case.As expected, the pre-collision profiles (see t = 1,800 s) do not yet differ substantially between the reference and channel configurations.At later times, the mechanical lifting appears more resilient in the wake of the merged front, owing to the frequent reflections which re-invigorate this region.This further illustrates how efficiently CPs can merge into stronger struc tures, as observed by Feng et al. (2015).Furthermore, mean and maximum CP heights are overall larger for the merged cases (up to 0.7 km at t = 14,400 s) and seem to decay very little over time.By contrast, the non-merged reference CP does not exceed 0.5 km height at any time and maximum heights appear to decay systematically over time.Thus, spontaneous triggering of new convective cells seems much more likely from merged CPs than individual ones.These observations remain valid when considering the 99th percentile instead.

Summary
This study investigates how key processes of CP currents, their propagation, collision and merging, depend on horizontal mesh resolution.By systematically refining the mesh while retaining equivalent CP forcing, we describe the evolution of the GF and the appearance of additional features and effects that depend mesh resolution.Our key findings are summarized below for each process.

Gust Front Propagation for a Single Cold Pool
The most basic configuration, simply referred to as propagation, focused on studying the formation and growth of a single CP in an idealized neutrally-stratified atmospheric boundary layer.Numerically, the flow was isolated from any outer perturbations and spurious reflections from the domain's boundaries.The forcing region, amplitude and time were chosen to be representative of a typical CP that might occur in an Radiative Convective Equilibrium (RCE) simulation.Therefore, the following results should be relevant to a wide range of atmospheric flow conditions.For a sensitivity study of the forcing zone size and amplitude, the reader is referred to Meyer and Haerter (2020).
The development of LC instabilities, typical of buoyancy-driven currents, was only observed at fine resolutions (Δx ≤ 100 m).Their onset was shifted further upstream as the resolution increased and occurred immediately after the end of the forcing period in the finest-resolution case (Δx = 25 m).To explore a possible causal effect between these two events, we ran a longer-forcing case to assess any correlation between them.It was found that the onset of the LC mode was virtually unchanged, though a continuous forcing strongly curtailed its growth rate-thus qualifying the nearly co-occurring events as coincidental.
Importantly, the LC mode acted as a converter of the efficient radial-vertical kinetic energy (i.e., contributing to front propagation and moisture updrafts) into the azimuthal mode.This fast transfer quickly depletes the leading GF's vortical strength, weakening the CP mechanical lifting.A visualization of the GF's cross-sectional rigid-rotational rate showed a noticeable drop for Δx ≥ 400 m.As the leading vortex flattened under the effect of the mesh coarsening, the amount of free-tropospheric moisture anomaly decreased.Free-tropospheric moisture is of particular interest as it determines the perturbation and moistening the troposphere right above the boundary layer upon passage of the front.If the moistening is strong enough trigger condensation, convection could be triggered.The GF traveling speed was found to be largely insensitive to the mesh resolution, although this specific configuration was constrained to a 10-km propagation radius.
Overall, a spatial resolution of Δx = 100 m appears sufficient to capture the LC-mode activation, mitigate numerical dissipation, and obtain a proper estimate of the mechanical lifting and GF depth.This is consistent with the findings of Grant and van den Heever (2016).Noteworthy, we observe that 99% of the 25 m-resolution upward moisture flux is resolved for Δx = 100 m, a value which quickly falters in coarser cases.The upward moisture fluxes exhibit an approximately self-similar profile across the [25-800] m-resolution range, following a power law of exponent −2.4,which could serve as the basis for developing a CP parameterization for weather forecasting purposes.

Gust Front Collision Between Two Cold Pools
This second configuration takes advantage of the domain boundaries' periodicity to split the forcing region into two-half cylinders, resulting into two fronts overlapping at the center of the domain.It is referred as the collision configuration.
We observed that a resolution of Δx = 400 m is sufficient to provide an accurate estimate of the total free-tropospheric moistening caused by the collision, both in amplitude and duration.Surprisingly, the peak in free-tropospheric water anomaly was shown to occur prior to the collision of the two vortex rings.Thus, new convective cells might sometimes be triggered ahead of gust fronts, rather than along with them.While it has been claimed that gravity waves cause such a phenomenon (Houze Jr, 2004), it is shown here to derive from the shape of the lower-tropospheric water anomaly bulge.Interestingly, the total water anomaly integrated above the boundary layer is, upon collision, almost twice as large as that caused by a single front propagation.Given that a collision comprises of two fronts-hence twice as much horizontal momentum-it is less efficient than a single front in converting horizontal momentum into vertical transport.This implies that some kinetic energy is lost to other modes during the collision.
While the distance-to-collision has been chosen to match the average CP-to-CP distance observed in RCE studies, the universality of these findings is unknown.Had the front traveled a longer distance before collision, the LC mode would have further developed and weakened the front prior to collision in the high-resolution case.
Conversely, an even shorter travel distance would have likely decreased the numerical dissipation observed for Δx = 800 m.

Merging of Multiple Cold Pools
The final configuration studied, referred to as merge, is a doubly-periodic channel with a longer forcing of 1h.
It enables us to study the coalescence and merging of a set of idealized CPs erupting simultaneously along the narrowest y-direction over a 4h-long runtime.The GF properties (speed and mechanical lifting) are quantified and compared to a reference collision-less case run on a (96 × 96 [km 2 ])-domain, which amounts to studying the same idealized CP as it propagates freely.This is the same 1h-long forcing case previously used as reference to evaluate the impact of a longer forcing time on the LC mode growth.This configuration was designed to quantify how merged fronts are able to strengthen each other, which was motivated by the recent findings from Jensen et al. (2022) showing how the simultaneous eruption and merging of CPs was instrumental to the formation of persistently-dry mesoscale regions in idealized simulations.
As in the propagation cases, we found little resolution dependence of the GF position after the first hour.However, by the end of the fourth hour, the final positions are ordered as a function of resolution, with the coarsest resolution lagging behind the finest one by nearly three km.This amounts to an averaged 5% slowdown in front traveling speed.We observed the finest Δx = [25, 50] m cases to be virtually indistinguishable.The comparison with the reference single-CP case highlights how the merged front is able to not only reach a higher traveling speed, but also sustain it for hours longer.This is a testament to the effectiveness of the CP merging process, as previously discussed (Feng et al., 2015;Jensen et al., 2022).Further, the merged front's upward moisture flux PDF was found to collapse remarkably well onto the same similarity profile, that is the −2.4 power law found for the propagation case.This robustness bodes well to numerical modelers as a unique flux-enhanced parameterization would benefit both CP processes, avoiding any complication in discriminating between them.Finally, we found that the flow in the wake of the merged GF is subject to further mechanical lifting due to reflected waves.Hence, convective activity could be better-sustained in this region as opposed to the wake of an isolated CP, which only leaves a subsidizing area of suppressed convection behind it.
To conclude, the merged front is stronger than an isolated CP, and it travels further and results in stronger and long-lasting perturbations of the free-tropospheric moisture field.It must be pointed out that these results are likely dependent upon the distance separating the CPs.Increasing (respectively decreasing) the distance between neighboring CPs would certainly diminish (respectively enhance) the merging between tangential gust fronts.

Closing Statements
The present study is a contribution to the expanding body of literature studying CPs, especially by numerical means.Using idealized large-eddy-simulations, we explore the effect of horizontal mesh resolution on the propagation, collision and merging of CPs.These three processes are of high-interest to weather and climate modeling as they enhance shallow circulation and are involved in the formation of MCSs.
Our findings suggest that a horizontal mesh spacing of Δx = 100 m effectively mitigates the excessive numerical dissipation occurring at coarser resolution and captures 99% of moisture upward motion.This provides a correct estimate of the free-tropospheric moistening by mechanical lifting, CP collision or incidental merging of several gust fronts.Δx = 100 m was also found sufficiently fine to accurately resolve the front's group velocity in any of the aforementioned processes.Interestingly, a robust −2.4 power law was derived.It could form the basis of a parameterization model aimed at compensating numerical dissipation by artificially upscaling the moisture fluxes in CPs.While we could not reach a numerical convergence on the growth rate and amplitude of the lobe-and-cleft instability, it appears clear that its development enhances the GF dissipation rate.Hence, further efforts are required to better understand how these effects would be included into a CP parameterization.
Although the mesh vertical resolution likely plays a crucial role in capturing the GF vortex, its impact on CP dynamics has not been considered in this study and would deserve further scrutiny.Likewise, much would be learned by extending the domain height up to the tropopause in order to capture the potential onset of CP-induced deep convection, and investigate its sensitivity to different ambient humidity rates.This would, however, inevitably increase the computational cost by another order of magnitude, strictly limiting the spatial resolution range of another numerical convergence analysis.Further, the CPs modeled in our study solely depend on a single triggering mechanism, both in size and amplitude.It is, therefore, unclear how our results would compare to more realistic and noisier configurations.Some elements of answers can be found in the previous study from Falk and van den Heever (2023) and Meyer and Haerter (2020).Finally, our results are certainly sensitive to the choice of numerical schemes.We expect higher-order, less dissipative schemes to perform better than the 5th-order WENO scheme in resolving the turbulent GF and the LC instability.

Figure 2 .
Figure 2. Time dependence of lowest level vertical velocity for all resolutions.Time and horizontal grid resolution are indicated along the horizontal and vertical axes, respectively.Colors ranging from dark blue to dark red indicate vertical wind speeds from −1 to +1 m/s.

Figure 1 .
Figure 1.Isosurfaces of potential temperature anomaly for single cold pool (CP).The isosurface corresponding to potential temperature anomaly θ′ = −0.1 K is colored by altitude z (km).The contours are ordered in time from left to right with times corresponding to images as indicated as ticks along the axis.The images offer a qualitative view of the idealized CP at various stages of its lifetime.See Movie S1 for an animation.

Figure 3 .
Figure 3. Appearance of lobe-and-cleft instabilities.Isolines of near-surface (z = 12.5 m) potential temperature anomaly θ′ = −0.1 K colored by local radial (outward) velocity u r .Each line corresponds to one timestep Δt = 30 s, such that the separation between two neighboring lines corresponds to the distance traveled during Δt.

Figure 4 .
Figure 4. Analysis of lobe-and-cleft instabilities for all resolutions.(a) Frontlines plotted in azimuthal space and colored by azimuthal velocity for Δx = [25, 50, 100, 200] m from top to bottom.The distance between each line corresponds to a 60 s increment (one-every-two timestep was discarded for the sake of clarity).(b) Single-sided amplitude spectra of the azimuthal velocity for all resolutions.(c) Measure of azimuthal asymmetry for all resolutions.(d) Mean frontline velocity for all resolutions.

Figure 5 .
Figure 5. Vertical displacement of water vapor for all resolutions.Azimuthally averaged contours of moisture perturbations plotted with velocity streamlines at t = 1,500 s for all Δx as indicated within the panels.Red crosses mark the locations of vertical velocity extrema.In each panel the dashed red line represents the vortex characteristic length which defines the front width.

Figure 6 .
Figure 6.Gust front (GF) width, water vapor perturbation and vortical strength versus resolution.Time-averaged (a) GF width δ and (b) time-averaged free-tropospheric (above z = 1 km) moisture gain G =  ∫  ′   as a function of horizontal mesh resolution Δx.The numbers correspond in (a) to the GF width in terms of the number of grid boxes (i.e., δ/Δx) and in (b) to the relative difference to the finest case (Δx = 25 m).(c) Profiles of leading vortex rotational rates plotted in time for all resolutions.The filled color markers correspond to the resolution case having the maximum vortical strength at a given timestep.Notice how the most-refined simulations are able to initially resolve a stronger vortex, which however dissipates at a faster rate than for the coarser simulations.Note that the horizontal axis direction is reversed (values decrease rightward) for all three panels.

Figure 7 .
Figure 7. Quantifying extremes in upward moisture flux versus resolution.(a) Probability density function of the positive upward moisture flux for all resolutions.(b-f) The same PDF(Δx) plotted versus their finer neighbor PDF(Δx/2).For each mesh, we determine the cut-off flux as the flux where the two lines depart from one another.This value is defined as the maximum flux bin where the condition ln(PDF(Δx)/PDF(2Δx)) < 0.25 is still fulfilled.(g) Cut-off fluxes as a function of resolution.Notice the inverse power-law trend.The numbers indicate the flux percentile resolved compared to the finest Δx = 25 m reference case.

Figure 8 .
Figure 8. Isocontours of potential temperature anomaly for gust front collisions.The isocontour corresponding to potential temperature anomaly θ′ = −0.1 K is colored by altitude z (see legend).The contours are ordered in time from left to right as marked along the axis and offer a qualitative view of the frontal collision between two idealized cold pool.See Movie S2 for an animation.

Figure 9 .
Figure 9. Water vapor anomalies under colliding gust fronts.(a) Cross-sectional contours of water mass fraction anomalies q′ in the y = 10 km plane at several times during collision for Δx = 25 m.(b) Integral of the z > 1 km water mass anomaly along the colliding (solid) and free (dashed) paths.(c) Maximum water mass anomaly in time for all resolutions.Note that the values along the free propagation are doubled for a more meaningful comparison.

Figure 11 .
Figure11.Comparison of gust front (GF) merging versus free propagation.(a) y-averaged GF x-location for the merged configurations and the reference 96 × 96 km 2 simulation run with Δx = 100 m where a single cold pool propagates freely for 4h (dashed cyan line).The front final positions are shown in the inset for all these cases.(b) The GF positions for the free-propagation and merging configurations (Δx = 100 m).Notice that the post-forcing growth rates resemble power laws and that the single front weaker speed drops from this power law earlier than the merged fronts.

Figure 10 .
Figure10.Isosurfaces of potential temperature anomaly for gust front merging.θ′ = −0.1 K colored by altitude z (km).The contours are ordered in time from left to right and offer a qualitative view of the cold pool as it merges with its neighbors, owing to the domain periodicity in the y direction.See Movie S3 for an animation.Note that only one half of the domain is shown.

Figure 12 .
Figure 12.Comparison of extremes in vertical moisture flux for merged and freely propagating cold pools (CPs).(a) Probability density function of the positive upward moisture flux in merged fronts for all resolutions.Notice the same scalability discussed for the single CP case in Figure 7.(b) Cut-off fluxes as a function of resolution.Notice a similarity profile with a different power law than for the single CP case.

Figure 13 .
Figure 13.Comparison of cold pool (CP) height profile for merged versus freely propagating CPs.Profiles of 90th-percentile front altitude, z 90 , plotted as a function of equivalent radius r* for the seven timesteps shown in Figure 10.The 100 m-resolution single-CP reference case is plotted in dashed line.Notice how the flow is less elevated in its wake (lower-altitude z 90 ).Using the 99th percentile instead yields similar results.