Detecting Cold Pool Family Trees in Convection Resolving Simulations

Recent observations and modeling increasingly reveal the key role of cold pools in organizing the convective cloud field. Several methods for detecting cold pools in simulations exist, but are usually based on buoyancy fields and fall short of reliably identifying the active gust front. The current cold pool (CP) detection and tracking algorithm (CoolDeTA), aims to identify cold pools and follow them in time, thereby distinguishing their active gust fronts and the “offspring” rain cells generated nearby. To accomplish these tasks, CoolDeTA utilizes a combination of thermodynamic and dynamical variables and examines the spatial and temporal relationships between cold pools and rain events. We demonstrate that CoolDeTA can reconstruct CP family trees. Using CoolDeTA we can contrast radiative convective equilibrium (RCE) and diurnal cycle CP dynamics, as well as cases with vertical wind shear and without. We show that the results obtained are consistent with a conceptual model where CP triggering of children rain cells follows a simple birth rate, proportional to a CP's gust front length. The proportionality factor depends on the ambient atmospheric stability and is lower for RCE, in line with marginal stability as traditionally ascribed to the moist adiabat. In the diurnal case, where ambient stability is lower, the birth rate thus becomes substantially higher, in line with periodic insolation forcing—resulting in essentially run‐away mesoscale excitations generated by a single parent rain cell and its CP.

Observationally, measuring cold pools is difficult, as high-resolution spatio-temporal records of the subcloud boundary layer are required to track a given spreading CP.Field campaigns thus often remain restricted to a particular, often spatially confined area (Feng et al., 2015;Hohenegger et al., 2023), use point-like or lowerdimensional measurements (Engerer et al., 2008;Zuidema et al., 2017), or need to resort to indirect signals, such as cleared areas in radar imagery (de Szoeke et al., 2017) or combinations of point and radar measurements (Kruse et al., 2022).Yet, from existing observational studies, it is evident that CP interactions help structure the cloud field and qualify as a mechanism to induce thunderstorm precipitation cells (Moncrieff & Liu, 1999).Measurements of the ocean temperature after CP events also point to lowered sea surface temperatures after pronounced CP occurrences (Pei et al., 2018), thus possibly giving rise to locally anti-correlated rain activity.
Numerical approaches hardly suffer from a lack of data coverage.Yet, in the past, three-dimensional simulations at cloud-resolving O(100m) resolution have not been possible over domain sizes of several hundred kilometers and weeks, required for reaching steady state dynamics.In recent years, however, the required scales have become accessible, and numerical simulations aimed at improving our understanding of cold dynamics have become feasible.For instance, process-oriented computer simulations have helped appreciate the role of cold pools in organizing the thermal and momentum structure within the subcloud layer (Böing et al., 2012).Further, numerical studies have helped elucidate the role of surface heat fluxes, induced by cold pools as they move along (Torri & Kuang, 2016), a dynamics which might affect the formation of so-called moisture rings (Chandra et al., 2018).The partitioning of surface heat fluxes into latent and sensible contributions was found to be strongly dependent on soil moisture, a variable which conversely also affects the size and propagation speed of cold pools due to changes in boundary layer humidity and rain shaft areas (Drager et al., 2020;Gentine et al., 2016).
On the larger-scale level of simulated mesoscale organization, cold pools have been suggested as a possible cause for non-random distributions of rain cells in space (Haerter et al., 2019), which may arise when CP gust fronts collide to trigger new cells at the loci of collision (Böing, 2016;Moseley et al., 2016;Torri & Kuang, 2019).Cold pools were also found to affect the moisture distribution in the subcloud layer (Böing et al., 2012;Schlemmer & Hohenegger, 2016) and may thereby aid the transitioning from shallow to deep convection (Kurowski et al., 2018).Even the paradigmatic "convective self-aggregation" (CSA), a system-scale symmetry breaking, was found to be affected by CP activity (Haerter, 2019;Jeevanjee & Romps, 2013;Muller et al., 2022;Nissen & Haerter, 2021).In an effort to study more realistic configurations, recent numerical experiments employ temporally varying surface temperatures to investigate the impact of the diurnal cycle on correlated CP activity, especially regarding the formation of mesoscale convective systems (MCS) (Haerter et al., 2020;Jensen et al., 2021), stimulating conceptual models to help explain the complex processes involved (Niehues et al., 2021).
Interestingly, only few attempts have been made at parameterizing cold pools in larger-scale climate models (Grandpeix & Lafore, 2010;Rio et al., 2009;Rooney et al., 2021).This may be due to an incomplete understanding of the fundamental processes affecting CP interaction, dynamics, and thermodynamic modifications of the boundary layer, as well as the pre-moistening of the lower free troposphere.Process-focused studies that leave out parts of the components involved (Fiévet et al., 2023;Haerter & Schlemmer, 2018;Meyer & Haerter, 2020;Romps & Jeevanjee, 2016), may help improve parameterizations further.
Several recent works have brought forward methods to track cold pools, in particular their gust fronts, in numerical simulations (Drager & van den Heever, 2017;Fournier & Haerter, 2019;Henneberg et al., 2020;Hirt et al., 2020;Rochetin et al., 2021;Schlemmer & Hohenegger, 2016).Yet, the detection of cold pools in numerical simulations is far from trivial: (a) the area affected by rainfall, where the CP is fed by rain evaporation, is often not cleanly delimited and subject to setting a threshold value, for example, on the precipitation rate; (b) the CP density current, as a highly dynamic object, experiences turbulent mixing and heat exchange with the surface and ambient environment, and is subject to spontaneous symmetry breaking, such as under the formation of "lobe-and-cleft" instabilities-again feeding back on the dynamics (Härtel et al., 2000;Markowski & Richardson, 2010;Simpson, 1972;Wakimoto, 2001); (c) the larger-scale pattern of cold pools and the rain cells produced by their collisions, is highly complex, with families formed by CP-rain cell networks making for a challenging tracking problem.
The current work addresses these challenges by building an algorithm which identifies cold pools and their causal chains based on a combination of buoyancy and dynamical measures.We apply our algorithm to cloud resolving simulation data based on radiative convective equilibrium (RCE) and diurnal cycle setups with and without wind shear to examine the impact of the simulation setup on CP dynamics.Based on the identified causal chains of cold pools, we investigate how the total accumulated rainfall of a given CP relates to (a) its capacity to trigger new rain cells and (b) the number of colliding cold pools which triggered its parent convection.Furthermore, we show examples of derived CP family trees and propose a simple conceptual model that replicates the observed cascadelike dynamics of cold pools within MCS.

Cloud Resolving Simulation Data
We employ the System for Atmospheric Modeling (SAM) (Khairoutdinov & Randall, 2006) to simulate the convective cloud field for a set of numerical experiments.The horizontal domain size is chosen to be square L × L, with L = 240 km in all simulations, and the horizontal grid spacing Δh is chosen as Δh = 200 m throughout, thus yielding N = 1,200 grid boxes in each horizontal dimension.Cyclic boundary conditions are chosen in both horizontal dimensions.We found this horizontal mesh resolution a suitable compromise, where relatively large domain sizes can be simulated for several days, yet, key CP effects, such as the dynamic gust front, can be simulated at satisfactory detail (Bryan et al., 2003;Fiévet et al., 2023;Hirt et al., 2020;Meyer & Haerter, 2020;Straka et al., 1993).The mesh is discretized along the vertical direction z using 64 levels of increasing depth Δz, ranging from Δz = 50 m at the first level (z = 25 m) to 1,000 m at the domain's top boundary (z = 26 km).The model resolves a non-hydrostatic anelastic form of the Euler equations, and uses the liquid water static energy and the total precipitating/nonprecipitating water mixing ratios as prognostic variables.We use the embedded onemoment microphysics scheme (Khairoutdinov & Randall, 2003) along with a first-order turbulent subgrid scheme (Smagorinsky, 1963).The radiative fluxes are calculated using the National Center for Atmospheric Research Community Atmosphere Model version 3 (Collins et al., 2006).
The surface conditions are prescribed and taken as a (saturated) sea surface of a given, horizontally-homogeneous surface temperature T s , that is, T s (x, y, t) = T s (t).We distinguish two sensitivities: 1. To mimic a diurnal cycle, we allow T s (t) to vary sinusoidally as T s (t) = T s ΔT cos (2π t/T).
(1) here, T s = 298 K and ΔT is either chosen as 0 K, which we term "rce," or as 2 and 4 K, respectively, termed "diu." 2. To simulate large scale wind shear, we additionally impose a height dependent wind tendency in the x-direction, with a piecewise linear profile given by (2) The velocity field is nudged toward this profile using a typical linear-relaxation term in the momentum equation with a timescale of τ = 1 day.This timescale is chosen to preserve the internal subdiurnal variability occurring through convective organization, such as through CP dynamics.
The above settings for ΔT and wind shear give rise to the six setups used in this paper, namely rce0K, diu2K, and diu4K, which can each have δ u = 0 or δ u = 16 m s 1 .For each setup we run the simulation for 7.5 days.
The integration time step, ∼ O(s), varies accordingly with a maximum admissible Courant-Friedrichs-Lewy (CFL) number of 0.8.The simulation output frequency is of one every Δt = 10 min.Each Δt, the following instantaneous fields are output: surface rainfall intensity I(x, y, t), near-surface specific humidity q(x, y, t) = q(x, y, z = 25 m, t), virtual temperature T v (x, y, t) = T v (x, y, z = 25 m, t), the horizontal wind vector u ≡ (u, v), with u = u(x, y, z = 25 m, t) and v = v(x, y, z = 25 m, t) as the lowest-level horizontal wind components in the x and y directions, and near-surface vertical wind speed w = w(x, y, z = 25 m, t).For simplicity, we refer to Δt as the "time step" throughout the paper.The discrete time t = t n ≡ nΔt, measured from the beginning of the simulation, is often denoted by the integer time step n.Similarly, we use integers i ∈ {1, …, N} and j ∈ {1, …, N} to label the two discrete positions x = x i ≡ iΔh and y = y j ≡ jΔh.Note that x i = x i+N and y j = y j+N in accordance with the cyclic boundary conditions.For a given 2D field Q we often use Q(i, j, n) to refer to the numerical value of Q at horizontal position (x i , y j ) and time t n .

Cold Pool Detection and Tracking Procedure
The main goal of the proposed algorithm is to detect and track cold pools in data from numerical simulations, while observing and identifying the involved causal relations between cold pools and rain cells: On the one hand, cold pools require evaporating precipitation and are thus inevitably linked to particular rain cells.To accommodate this causal relation, rain cells are taken as the origin of cold pools in the proposed algorithm.Thus, a rain cell is a required condition for the identification of a CP.Each rain cell in the domain is assigned a unique rain identification number (ID).In case a rain cell generates a CP, the CP receives the same ID as the associated rain cell.On the other hand, CP gust fronts can trigger new rain cells, which might again generate new cold pools (Böing, 2016;Haerter et al., 2019;Tompkins, 2001a).Accordingly, the algorithm stores rain cells triggered at a CP gust front, as well as potential cold pools emerging from these rain cells, as children of that particular CP.In this manner, the algorithm not only detects and tracks cold pools, but also their relationships.
The proposed CP detection and tracking algorithm (CoolDeTA) consists of four main steps: (a) identification of rain patches (RPs), (b) segmentation of the domain into "potential CP" and "no CP" regions, (c) derivation of markers, and (d) labeling of cold pools.These steps are performed for every output time step of the respective simulation or simulation time window.
For notational convenience we introduce several definitions: • The deviation of any given scalar quantity Q(i, j, n) relative to its spatial mean is denoted as where Q(n) denotes the spatial mean at t n .• Occasionally, we require Q at a given time step t n to be normalized to the range between zero and one.This is accomplished by the rescaling where min(Q(n)) or max(Q(n)) refer to the minimum and maximum of the two-dimensional field Q(n) at time step t n .• For simpler notation we write the Kronecker delta symbol as δ(a, b), which equals unity for a = b and zero otherwise.• Two (spatial) indices (i, j) and (i′, j′) are considered contiguous, if dist((i,j) (i′,j′)) = 1, where dist denotes the shortest distance operator, which ensures cyclic boundary conditions on the toroidal square lattice geometry.• We use p r (i, j, n) and p cp (i, j, n) to indicate N × N integer fields for labeling RPs and CP patches at time step t n .
Analogously, we also define an auxiliary field p seg (i, j, n).For either, we define a superscript m as labeling the subset of indices (i, j, n) of p r or p cp , where the integer is equal to m, for example, (i) Identifying RPs.At any discrete time t n a RP is defined as a spatially contiguous region in which the rain intensity I(x i , y j , t n ) ≥ I 0 .Generally, we employ I 0 = 2 mm hr 1 , which has been determined as critical rate for downdrafts penetrating into the subcloud layer (Barnes & Garstang, 1982).The algorithm assigns a unique ID to every RP identified by this procedure.In this way, p r (i, j, n) is populated at each time step t n , where p r (i,j,n) = { 0 if I(i,j,n) < I 0 otherwise: ID of respective rain patch. (5) where s′(i, j, n) are the horizontal wind speed fluctuations, computed as I seg (i,j,n) is a measure which we assume to be increased in areas where cold pools are present.Since cold pools are associated with fluctuating horizontal wind speeds but decreased virtual potential temperatures, T′ v (i,j,n) in Equation 6is subtracted from unity, that is, relatively low T v increases the probability of a CP.The two quantities are combined with equal weights to ensure equal importance is given to both thermodynamic and dynamical features.We choose perturbations rather than absolute quantities in Equation 6 to be able to apply the method to simulations with other boundary conditions, such as non-zero surface winds, and thus increase the method's scope of application.
A k-means algorithm (Pedregosa et al., 2011) now subdivides every grid cell (i, j) into two cluster groups, namely "potential CP" or "no CP" by minimizing the one-dimensional distance metric within both cluster groups.The result is a horizontal N × N binary segmentation p bin (i, j, n) with p bin (i, j, n) = 1 for all "potential CP" grid cells and p bin (i, j, n) = 0 for all "no CP" cells.As an example, Figure 1a shows the corresponding histogram of I seg for diu4K at time step t 471 representing a morning scene.The lower and upper cluster groups represent regions where cold pools are considered prohibited versus feasible.A "cold pool" grid cell must fall within the feasible cluster group.Yet, additional criteria below are required.
Spreading quickly beyond the boundaries of their RPs, fully developed cold pools should at least encompass the same area as their related RPs.Thus, we apply the minimum area A 0 criterion from (i) also to the identified "potential CP" regions.For this reason, spatially contiguous "potential CP" regions in p bin (i, j, n) are identified and labeled with a unique integer number.The result is, analogous to (i), a horizontal N × N integer field p seg (i, j, n) where a patch, labeled m, is only retained if its area (iii) Deriving markers.A marker could be defined as a center of mass (COM) of a given RP m at time t n .We however found that the actual (dynamical) downdraft within this patch more precisely marks the CP center.A good choice for the downdraft field I Tw (i,j,n) was found to be a combination of normalized temperature and vertical velocity, namely so that for a given patch, labeled m, which is large for locations of minimal temperature and vertical velocity.The corresponding COM, (C m x ,C m y ) , follows as  (i, j, n), where p mar i m mar ,j m mar ,n) = m and which is zero for all locations without markers.
When the rain associated with a CP labeled m from the previous time step t n 1 stops, the CP will have no active RP p m r in p r (i, j, n) anymore, that is, ∑ N i,j=1 δ p r (i,j,n),m) = 0.However, cold pools may still spread and trigger new convection.Therefore, we in this case introduce their latest rain marker to p mar (i, j, n) instead.Apart from the rain marker in the center of the active or latest RP, cold pools which were present in the previous time step t n 1 obtain an origin marker.An origin marker represents the first rain marker which initially generated the CP, that is, the center of its initial RP.Analogous to the rain markers, origin markers are introduced to p mar (i, j, n).The combination of origin and rain markers increases the robustness of the final CP labeling, which will be illustrated in step (iv).
(iv) Labeling cold pools.For a time step t n , the final field with labeled cold pools, p cp (i, j, n), is-analogous to the RP field p r (i, j, n) (compare Equation 5)-a horizontal N × N integer ID field where p cp (i, j, n) specifies the ID of a CP or p cp (i, j, n) = 0, if no CP is present at (i, j) and time index n.The CP field is obtained using a watershed algorithm (van der Walt et al., 2014), which can be thought of as placing different water sources in a topographical map and flooding it up to the highest elevation (see Figure 1c), such that the water table is equal in all basins.The resulting watersheds represent the boundaries between the different water sources and the corresponding basins represent regional minima.
In order to derive the labeled CP field, we provide the watershed algorithm with three ingredients: the locations of the water sources, that is, the marker locations p mar (i, j, n) computed in step (iii), an elevation map I top , and a mask, which acts like a stencil preventing certain regions from being flooded.The elevation map is computed as Cold pool regions are characterized by cold (low T(i,j)) and often dry air (low q(i,j)) which forms basins with regional minima in the elevation map, whereas horizontal convergence with associated vertical updrafts towers up walls at CP boundaries (large ŵ(i,j)), thus separating the CP basins from each other.The mask, as a final ingredient to the watershed algorithm, is provided by the binary segmentation p bin (i, j, n) from step (ii).Whereas CP1 and CP3 in Figure 1c are separated by the mask's "no CP" region between them, CP2 and CP3 are in contact with each other.The contact location, that is, the point where the two cold pools collide (in 2D, it is a line), depends on the elevation map I top .Since water from different sources, representing different cold pools, does not mix, the contact location remains constant for the rest of the flooding.
Since only its gust front should separate a CP from its environment, we do not allow holes, that is, enclosed grid cells with different ID values, within individual cold pools.Thus, we close potential holes within individual cold pools by assigning each grid cell of the hole the ID of the surrounding CP.By this procedure, we thus merge cold pools if the grid cells of the hole carry the ID of other cold pools.
In later stages, a CP basin in the elevation map can consist of multiple local minima, caused by newly formed RPs within the CP area or an advected initial RP.By providing both rain and origin markers to cold pools which were present in the previous time step t n 1 , we make sure that cold pools in later stages of their life cycle are not accidently flooded by neighboring cold pools.
Figure 1b summarizes the main steps of CoolDeTA.Note that both tracking and the interactions between steps (ii) and (iii) will be described separately in the following Sections 2.3 and 2.6, respectively.

Cold Pool Tracking
Provided that RPs, p r (i, j, n 1), and/or cold pools, p cp (i, j, n 1), from the previous time step t n 1 are available, we track cold pools already during the identification of RPs in step (i).Since a RP, labeled m in p r (i, j, n), obtains a corresponding marker p mar i m mar ,j m mar ,n) = m, an associated CP which might be identified based on the watershed algorithm in step (iv) would be labeled m in p cp (i, j, n), as well.The label m thus represents the ID of a RP in p r (i, j, n) and the ID of an associated CP in p cp (i, j, n).Accordingly, cold pools can be tracked by adjusting the label of tracked RPs in step (i).For this purpose, we evaluate the spatial overlap of patches from consecutive time steps.A RP labeled m′ at time step t n is tracked if the overlap with a RP labeled m from time step Journal of Advances in Modeling Earth Systems where 1 is the indicator function and the set where the set , the previous CP provides its label, that is, the label of RP m′ is replaced with the label m 2 .In case RP m′ is overlapped by multiple RPs from the previous time step with deviating labels, the RP with the largest overlap, O max r , provides its label.For this reason, the new RP in time step t n of Figure 1d which overlaps RP1 and RP2 from the previous time step t n 1 is tracked as RP1.
Apart from the rain markers derived from RPs, all cold pools in p cp (i, j, n 1), which existed in the previous time step, receive origin markers and-if they do not have an active RP anymore-also their latest rain marker (see step (iv)).Accordingly, CP1, CP3 and CP4 in Figure 1d each receive two markers in time step t n , representing the COM of the current RP and that of the initial RP from the previous time step t n 1 , respectively, whereas CP2 only receives its latest rain marker from time step t n .However, due to its origin marker CP2 is still present at time step t n although it does not possess an active RP anymore.

Cold Pool Relationships
The association of RPs and cold pools allows the method to consistently label cold pools.Beyond this, through the spatial dependencies, we are also able to capture relationships between individual cold pools.By tracing which CP (gust front) triggered which RP and thus also a potential CP resulting from it, CoolDeTA reveals CP networks and their associated causal chains.
In order to describe the relationships between cold pools, we introduce the roles of "parent" and "child."A RP labeled m′ at time step t n which was not present in time step t n 1 , is a child of a CP, labeled m, from time step , that is, if the RP intersects the boundary of the CP, but is not fully located within its area.Conversely, CP m becomes a parent of RP m′.If RP m′ is a child of CP m and generates a new CP m′, the new CP adopts CP m as a parent and becomes its child, as well.Since cold pools can trigger new rain events, they can become parents-both of RPs and the resulting cold pools-as well as children, whereas RPs can only be children.
In case a RP labeled m′ at time step t n which was not present in time step t n 1 overlaps multiple cold pools from time step t n 1 , any CP m with an overlap 0 < O cp (m,m′,n 1,n) < f cp × A p m′ r ) becomes a parent of RP m′ and thus also of potential cold pools resulting from it.

Cold Pool Dissipation
When the rain associated with a CP stops and thus no longer feeds the density current anymore, the surface boundary layer recovers until the CP is eventually dissipated.However, providing origin markers to cold pools from the previous time step t n 1 in time step t n , cold pools could potentially exist forever, even after their rain has stopped.To rule out infinite lifetime, we implemented a dissipation mechanism within step (iii) of CoolDeTA.To this end, we differentiate between "dissipating" cold pools that are still in the process of dissipating and "dissipated" cold pools that have completely dissipated.We define a CP, labeled m, as dissipating at time step t n if it has no active rain anymore, that is, ∑ where the set For each CP we count the number of time steps in which it is dissipating, n dis .We consider a CP, labeled m, as dissipated if n m dis ≥ 3 and O seg (m, 1, n last , n) = 0, that is, if it is dissipating for ≥3 time steps (=30 min) and if its latest RP p m r is completely classified as "no CP" in p bin (i, j, n).A CP which is dissipated at time step t n does not obtain markers anymore in p mar (i, j, n) and will thus not exist in the labeled CP field p cp (i, j, n).

Additional Algorithm Rules
The evolution of the populated fields of both cold pools, p cp , and RPs, p r , is highly complex.On the one hand, relatively cool air associated with remnants of dissipated cold pools might lead to the misclassification of grid cells as "potential CP" in the segmentation, p bin , in time steps with very weak or suppressed convection and thus only weak down and updrafts.On the other hand, processes such as the merging of RPs or cold pools, as well as the formation of new RPs at CP gust fronts, complicate the tracking.In the following, we introduce several rules implemented in CoolDeTA to increase its robustness with respect to these and other cases.

Divergence Criterion
Cold pools are driven by atmospheric density gradients resulting from the evaporation of rain.Evaporative cooling causes the affected air to sink toward the ground and spread outwards, creating wind gusts.Consequently, cold pools are associated with a horizontally divergent flow in their center and a horizontally converging flow at their gust front.Thus, we only keep a "potential CP" region, labeled m in p seg (i, j, n) at time step t n , in the segmentation, p bin (i, j, n), if the mean divergence of the horizontal wind field, ∇ ⋅ u(i,j,n) ≥ 0 for grid cells (i, j) in its patch interior int p m seg and if ∇ ⋅ u(i,j,n) ≤ 0 for grid cells in its patch boundary ∂p m seg .The divergence is approximated using a central difference scheme.We define a grid cell (i, j) as part of the interior int p m seg of patch m if all contiguous grid cells are contained in p m seg .A grid cell (i, j) is part of the boundary ∂p m seg of patch m if it encompasses at least one contiguous grid cell contained in p m seg and one contiguous grid cell not contained in p m seg .Since the "potential CP" region m in p seg (i, j, n) can contain multiple CP instances and thus colliding gust fronts, we exclude a grid cell (i, j) where σ d represents the standard deviation of the horizontal wind field divergence and z 95th = 1.645.
Furthermore, new rain cells might be triggered at the boundary of the "potential CP" region m compensating its convergent flow.Analogous to int p m seg we thus exclude a grid cell (i, j) ,j,n) ≥ 1 mm/h with z 75th = 0.675.We chose deviating factors z 95th and z 75th for the exclusion of grid cells from CP interior and boundary, as colliding CP gust fronts will be associated with relatively strong convergence in the patch interior, whereas new rain cells triggered at the boundary will be associated with relatively weak divergence compared to the domain mean.The values of both z 95th and z 75th were determined by experiment.The exclusion of potentially misleading grid cells, as well as the weak conditions in terms of int p m seg not being convergent and ∂p m seg not being divergent ensure that no correctly classified "potential CP" region is dropped.
Due to CP interactions such as collisions and the formation of new cold pools at gust fronts, a CP can loose its occupied area to another CP.To not supersede this process artificially, a CP, labeled m, at the previous time step, t n 1 , only obtains an origin marker p m mar i m mar ,j m mar ,n) = m if it still occupied the marker location i m mar ,j m mar ) at t n 1 , that is, if p cp i m mar ,j m mar ,n 1) = m.
Apart from CP interactions, where a CP gust front can be pushed toward the CP center, CP gust fronts spread outwards only, due to the horizontally diverging flow in the CP interior.As long as a CP, labeled m, is not dissipated and receives at least one marker at time step t n , we thus adapt the segmentation p bin (i, j, n) in step (iii) so that all grid cells occupied by the CP at the previous time step t n 1 , p cp (i, j, n 1) = m, are labeled as "potential CP" in p bin (i, j, n).Whether CP m looses parts of its occupied area to an intersecting CP in time step n depends on the topography of the elevation map, I top (i,j,n), and is determined during the watershed flooding in step (iv).

Additional Overlap Criteria
When new cold pools start to form, grid cells with strong negative temperature perturbations inside the downdraft region might already be classified as "potential CP" in p bin .In order to detect fully developed cold pools only, a new RP, labeled m, detected at time step n, obtains a rain marker at time step The final check is applied to the labeled CP field p cp derived based on the watershed flooding in step (iv).By comparing the area, A cp , of newly formed cold pools with the area, A r , of their RPs, we get an indication of erroneous labeling.Erroneous labeling can result from cool air remnants being classified as "potential CP" during time steps with very weak or suppressed convection.In this case, markers of isolated rain events could flood areas A cp ≫ A r .Another potential source of erroneous labeling is new RPs, triggered at CP gust fronts: Since not each of these RPs generates an independent child CP which is separated from its parent CP by their colliding gust fronts, child cold pools may flood into the area of their parent cold pools and vice versa.Thus, an independent and fully developed CP should, on the one hand, at least feature the size of its RP.On the other, a new CP cannot occupy an area substantially larger than its RP, given the time step of Δt = 10 min.To intervene only in erroneous cases, we keep a new CP, labeled m, if 1 ≤ A(p m cp )/A p m r ) ≤ 3. Otherwise, we either drop it, that is, p m cp = 0, if it has no other CP in its neighborhood, or we assign the label of another CP, labeled m′, that is, p m cp = m′, otherwise.In the latter case, the CP label m′ is determined based on its number of parents, k p : if k p = 0, it is either the label of the CP which overlaps the largest area of RP m (if A(p m cp )/A p m r ) < 1), or the label of the CP which lost the largest area with respect to the previous time step CP m is simply assigned the parent label, m p , that is, m′ = m p .The final case, k p > 1, is analogous to the first but based on the parents of CP m.
Table 1 summarizes the parameters and conditions used in the proposed algorithm.Recommendations for the application of CoolDeTA to other simulation setups are provided in Text S1 in Supporting Information S1.

Analysis of the Simulation Data
We consider the first 3 days (432 time steps) of the simulations as spin-up period.The spin-up period is chosen so that the total domain rainfall per time step, R, has visibly reached a steady state in both the RCE and diurnal cycle setups (Figures 2e and 2f).Further time series are visualized in Figure A1.To be able to track the full life cycle of cold pools and their causal relationships, the onset of convection is an ideal starting point for the analysis with CoolDeTA.Apart from rce0K, where R is essentially constant (Figure 2f), we identify the corresponding starting time step based on the standard deviation of the virtual temperature, σ Tv , and its time derivative, σTv .Based on their 1 hr running average, we find σ Tv ≥ 0.15 K and σTv ≥ 0.02 K s 1 as good indicators for the onset of convection in our simulations.In the diurnal cycle setups, as well as the rce0K setup with wind shear, we choose the starting time step on day four, t start , so that the two conditions are fulfilled at least for the next hour (6 time steps).For rce0K we simply start the analysis with CoolDeTA at the first time step of day four, at t 432 .
In the diurnal cycle setups, R oscillates nearly periodically (Figure 2f).Especially in the two diu4K simulation setups, convective activity is largely suppressed during time steps with low T s .Since the k-means algorithm always clusters the horizontal domain into two groups, namely "potential CP" and "no CP," the segmentation cannot work properly in the absence of cold pools.In such cases, also the additional rules may be insufficient.For this purpose, we stop the CoolDeTA analysis if R px < 0.001 mm for more than 30 min (3 time steps) or as soon as σ Tv < 0.15 K ∧ R px < 0.0015 mm, with R px = R N×N .Analogous to the definition of the start time step, we apply the conditions to the 1 hr running averages of σ Tv and R px .The stop conditions affect only the two diu4K setups, which we stop at the end of each day and re-start on the following day based on the defined start conditions.

Results
In the following, we present results of our CP detection and tracking method.We first compare our method to previous ones (Feng et al., 2015;Gentine et al., 2016;Tompkins, 2001a).We then employ our current method to discuss the effect of the simulation setup, that is, wind shear or diurnal cycles, on CP expansion.Third, we make use of our method to analyze the causal relationships between cold pools and their effect on CP rainfall.Based on our findings we then propose a simple model to capture CP spreading-triggering dynamics.We finish the chapter by visualizing the family trees of two cold pools.

Comparison to Other Methods
Several previous CP detection methods have exploited thermodynamic CP features (Feng et al., 2015;Gentine et al., 2016;Tompkins, 2001a).While Feng et al. (2015) and Tompkins (2001a) use buoyancy thresholds of 0.005 m 2 s 1 and 0.003 m 2 s 1 , respectively, for the detection, Gentine et al. ( 2016) apply a k-means algorithm to the horizontal virtual temperature field of the first model level to distinguish cold pools from their environment.For an RCE-like setup, rce0K (Figure 2a), as well as a diurnal cycle case, diu4K (Figures 2b-2d), we compare CoolDeTA to such existing methods.The time steps shown are highlighted in the time series (Figures 2e and 2f) as thin vertical lines.Inspecting the different examples (Figures 2a-2d), it is apparent that the CP dynamics in rce0K occurs at relatively small scales O(10 km), whereas the dynamics in diu4K ranges from scattered small scale events during the early stages of the diurnal cycle to elaborate mesoscale features O(100 km) at the late stages of the diurnal cycle.
A primary goal of CoolDeTA is to detect both thermodynamic and dynamical features of cold pools.To obtain a reference for the latter, we highlight grid cells of pronounced positive vertical velocity in each panel, typically indicating locations of dynamical CP gust fronts (red pixels in Figures 2a-2d).Inspecting the plots, it is apparent that on several occasions, the detection results by all four methods nearly coincide, especially during early stages of CP expansion.However, important differences can also be seen, notably within the inset to panel (a), where CoolDeTA is able to capture the locations of the dynamic gust front (red pixels) reliably, whereas methods based on thermodynamics typically indicate gust fronts that are located further inward, seen for example, by comparing the lines of different colors in the inset.As is visible on the left side of the larger CP patch in t 473 of the inset, CoolDeTA keeps dissipating CP instances as long as they feature a partly active dynamic gust front.
For the early stage of the diurnal cycle (Figure 2b) the agreement of CoolDeTA, Gentine et al. (2016) and Tompkins (2001a) is arguably good.However, as cold pools spread in the course of the day, the dynamic gust fronts tend to detach from the thermodynamic ones, leading to an increased discrepancy between the methods.It can be visually verified that CoolDeTA is still capable of identifying most regions of strong vertical velocity.
Apart from the common minimum area A 0 which we apply in Figures 2a-2d for CP patches of all methods, CoolDeTA applies the threshold A 0 also to RPs associated with individual cold pools and conducts further checks to confirm the existence of a CP (see Section 2.6).Consequently, the minimum CP area detected by CoolDeTA tends to be larger than for other methods.As a case in point, some of the smaller patches identified by the other methods either suddenly disappear or do not feature a spreading dynamic gust front (Figure 2a, enlarged area), suggesting that CoolDeTA manages to successfully detect only robust cold pools.

Effect of the Simulation Setup on Cold Pool Expansion
We now employ CoolDeTA to investigate CP temporal expansion and how it depends on the choice of simulation setup.Cold pool expansion is a basic dynamical CP feature, which we here capture as the change in CP area, A cp , over time, on the one hand for cold pools without offspring during their lifetime (Figure 3a) and on the other for cold pools with offspring (Figure 3b).To not distort the results with cold pools which shrink due to CP interactions, we only consider a CP as long as its area is expanding and thus A cp (n) > A cp (n 1).
Conditioning on cold pools without children (Figure 3a), it is noticeable that within the log-log representation, all curves start out following a modest increase of CP area with time.To make contact with previous works suggesting power law dependencies of CP radius on time (Meyer & Haerter, 2020;Rooney, 2015), we indicate a dependency of A cp ∼ t α with α = 1.2 (dashed gray lines).After approximately 1 hr of expansion, the curves of setups without wind shear somewhat increase in exponent, thus departing from a power-law dependency.We attribute this to the higher density of cold pools for setups without wind shear, which increases the probability for cold pools to merge with others.The departure is even stronger for diu2K and diu4K where cold pools increasingly cluster within thermodynamically favorable subregions.
The curves for cold pools with children exhibit an increased power-law exponent already at earlier stages (Figure 3b).In comparison to the curves of cold pools without children, the exponent seems to be increasing more strongly, especially during the later stages of expansion.We interpret this strong increase as a consequence of enhanced triggering near CP edges for cold pools with children, whereby the original CP can merge with additional offspring excited by it-thus developing into a mesoscale CP.Apart from this, the diurnal cycle again appears to promote rapid CP expansion, especially during the later stages.
In both Figures 3a and 3b, cold pools under wind shear conditions have larger areas from the outset.A possible reason is the overlap that a RP is required to have with a "potential CP" region in order to obtain a marker and thus form an associated CP: when wind shear is present, rain might reach beyond the emerging CP and limit their overlap.

Effect of Cold Pool Relationships on Cold Pool Rainfall
Given these findings on CP dynamics we now investigate: (a) how the total accumulated rainfall corresponding to any given CP relates to its capacity to trigger offspring (Figure 4a); and (b) how this rainfall relates to the number of parents of a given CP (Figure 4b).For this purpose, we define the total accumulated rainfall of a CP labeled m as Overall, the majority of cold pools do not give rise to offspring, with approximately 50% of cold pools without detectable children (Figure 5a).However, as we inspect the probability distributions, it is found that, moving from rce0K to diu4K, the distribution function is close to an exponential for rce0K (dotted line), whereas it is well approximated by a power-law ∝ k 1.5 in the case of diu4K (dashed line).This qualitative finding is also mirrored by the wind shear cases.A power-law distribution hints at a "rich-gets-richer" feedback, where cold pools that have already grown, may be more likely to grow further.
At the same time, it is clear that for large numbers of CP children, associated RPs tend to have much greater total rainfall rates (Figure 4a).These rates range from a factor of two for rce0K to more than an order of magnitude for diu4K.Notable differences also exist between the different setups: diurnal T s -forcing and wind shear seem to promote larger numbers of children and generally come with greater parent rainfall.
Apart from cold pools with one parent, the relation between parental rainfall and the number of children is almost perfectly mirrored by the number of parents (Figure 4b): larger total rainfall typically implies a larger number of parents.Also, for setups with T s -diurnal cycles or wind shear, the role of parent cold pools contributing to rain cell formation increases.Yet, the rainfall rates of cold pools with one parent hardly differ from those of cold pools without a parent.We attribute this to the higher risk that a single-parent CP has of getting merged into the parent CP: whereas multiple parents are associated with CP collisions, interrupt the expansion of the parents, single-parents might still spread and thus merge with excited offspring.

Simple Cold Pool Offspring Model
To capture CP spreading-triggering dynamics, we suggest a simple model (Figure 5b) where cold pools primarily grow by peripheral spawning of new rain cells, thus child cold pools, that then help the CP expand further.The model is hence different from those in the recent literature (Böing, 2016;Haerter, 2019;Haerter et al., 2019;Nissen & Haerter, 2021) in that the dynamics described previously was inward-directed, such that existing cold pools could excite new convection by collisions among two or three gust fronts.As a result, a given population of cold pools could topologically only cause new convection within the interior of an ensemble of spreading cold pools.In the current model we describe an outward-directed dynamics, where a given evolving mesoscale convective system, embedded in an unstable environment, can give rise to cascading additional convective activity at its periphery-thus leading to an outward expansion of the population of convective cells.
In our model each CP is initialized with an area a of a = a 0 ≡ 6 km 2 representing the median area of all cold pools from Figure 3 at t cp = 10 min, and circumference l = 2(πa 0 ) 1/ 2 ≈ 8.7 km.This circumference is subdivided into segments, each of length l 0 ≡ 1.6 km, taken to represent the approximate diameter of a new rain cell of area A 0 .
Within each segment and time interval Δτ ∼ 20 min, chosen to represent the timescale of a convective rain cell, the model allows for a new rain cell to be spawned with a small probability 0 < p ≪ 1.The constant l 0 ≪ l is thereby assumed to be substantially less than the total circumference l, such that even small cold pools can have several segments.Any new rain cell will add an area increment a 0 to the existing CP area a.In addition, as the CP area increases, its total number of segments n s = l/l 0 will also increase due to the larger circumference.Thus, larger cold pools will typically have a larger number of successful spawning events.
As mentioned, the dynamics proceeds in discrete time steps Δτ, during which each segment will be able to experience spawning at probability p.If none of the segments give rise to spawning during the time step, the CP is considered terminated and no further expansion is possible.If a number m s > 0 of segments give rise to spawning, the CP area is incremented by Δa ≡ m s a 0 and the expansion proceeds iteratively (tree diagram in Figure 5b).At each fork of the tree diagram, at probability 1 P(a i ) the CP expansion is terminated. 1 P(a i ) represents the probability that none of the segments spawn a new rain cell, that is, for n s = l/l 0 segments and spawning probability p, 1 P(a i ) = (1 p) n s . (16) Conversely, the probability for continued CP expansion is As cold pools expand, n s continues to increase and P(a i ) → 1-thus allowing larger cold pools to experience a larger probability of further expansion.
Simulating the model using a total of 30 time steps Δτ and 10 5 CP realizations allows us to mimic the statistics found for the cloud-resolving simulations (Figure 5c): whereas small values of spawning probability, for example, p = 0.035, yield a nearly exponential distribution of the number of children, a larger value of p = 0.08 gives rise to approximate power-law decay.Notably, as cold pools grow to large areas, spawning will be essentially guaranteed at subsequent time steps-thus giving rise to run-away expansion.This is verified by allowing for a greater number of 60 time steps (thin brown curve), where a fraction of "successful" cold pools accumulates ever larger numbers of children.Using this model, our interpretation of the CP dynamics is twofold: spawning probability p may be determined by the convective instability surrounding a given parent CP, and different p can yield very different offspring distributions.
For the RCE case, p may be overall modest, as the atmosphere tends to be in a marginally stable state close to the moist adiabat.For the diu cases the situation is quite different, where activity is invigorated by the destabilizing effect of the increased diurnal surface heating-thus allowing for larger values of p, yielding potentially longtailed offspring distributions.A complete run-away effect is likely prevented by the finite duration of the diurnal heating.It may thus be interesting to revisit the statistics (Figure 5a) for an artificially lengthened diurnal cycle, say, corresponding to a 48 hr-day.

Cold Pool Family Trees
Based on a subsection of the full simulation domain we now visualize an example of the temporal evolution of a labeled CP field from diu4K (Figures 6a-6c), along with the corresponding family trees of the two initial cold pools (Figure 6d).The subsection of interest is indicated as black square in the inset which shows the specific humidity anomalies, q′, for the full simulation domain at t 764 .Like most of the other regions in which cold pools start forming (dark red blobs with large gradients), the subsection of interest features a moisture-rich surface layer.Here, the two initial cold pools depicted in Figure 6a, CP46 (yellow patch) and CP49 (blue patch), formed at t 763 and t 761 , respectively (Figure 6d).
In Figure 6b, 80 min later, the two cold pools have expanded and several offspring have been generated.However, CP848 (red ID in Figure 6d), one of the first children of CP49 which was triggered at the freely expanding gust front, could only exist for one time step until it was merged by its parent.Although CP848 contributed to the formation of the new CP CP1169 (see Figure 6d), it is thus not present anymore at t 772 .Comparing Figures 6b and  6c, it is apparent that also the other cold pools which were triggered at freely expanding gust fronts, namely CP684, CP982, CP1169 and CP1172, shared the same fate.Only the two cold pools CP466 and CP835 in the collision zone of CP46 and CP49 are still present at t 780 and managed to expand over the former area of CP46 and CP49.Since it is not evident from Figure 6b we note that CP466 and CP835 formed before the collision of CP46 and CP49, otherwise both CP46 and CP49 would be their parents.
A number of CP detection methods exist, which often make use of a buoyancy threshold (Feng et al., 2015;Gentine et al., 2016;Tompkins, 2001a) or a detection of the gust front through a criterion on the dynamics (Fournier & Haerter, 2019;Henneberg et al., 2020).We have presented a method, CoolDeTA, that combines thermodynamic and dynamic CP signatures.CoolDeTA detects and tracks cold pools as well as the location of their dynamical gust front.On this basis, CoolDeTA is able to identify relationships between cold pools and subsequent rain cells and thus CP "offspring."CoolDeTA therefore can be used to analyze interactions between one or more cold pools and the ones emerging in the surroundings.
Like most other methods (Drager & van den Heever, 2017;Feng et al., 2015;Fournier & Haerter, 2019;Gentine et al., 2016;Henneberg et al., 2020), CoolDeTA cannot fully refrain from setting some threshold values.In the present work we apply a surface rainfall intensity threshold, I 0 = 2 mm hr 1 , to delineate spatially contiguous RPs, and a corresponding area threshold, A 0 = 2 km 2 , for RPs and potential CP regions.Thresholds are generally case dependent and thus limit a method's scope of application.Although particularly the CP labeling of RCE frameworks would occasionally benefit from an even smaller A 0 , we however claim that both of our thresholds serve only as a preselection of robust rain and CP patches and do not impact the final CP contours: Since we provide each new RP above A 0 one marker in its downdraft COM, all of these RPs have the same chance of forming a CP.The CP contour, and whether a RP forms a CP at all, depends on the result of the k-means algorithm.A sufficiently large RP and an associated marker are accordingly necessary, rather than sufficient, conditions with respect to the formation of a CP.Also the additional conditions summarized in Table 1 serve solely the purpose of optimizing the preselection in order to increase the robustness of the CP labeling.To verify this, we deactivated the dissipation threshold, n dis , as well as all additional conditions and reproduced Figures 2 and 6 based on this parameter setup in Figures S1 and S2 in Supporting Information S1, respectively: whereas the deactivation causes essentially no change in the diu4K plots (Figures S1b-S1d in Supporting Information S1), several of the smaller patches in the rce0K plot which were identified as cold pools by other methods but not by CoolDeTA before, get classified as cold pools without the additional conditions (Figure S1a in Supporting Information S1).The same applies to the CP relationships (Figure S2 in Supporting Information S1), which are identical under deactivated conditions, except for a few additional CPs that disappear shortly after being detected.With respect to dissipating cold pools, the omission of n dis can even positively affect the results in the rce0K plot since some of the weak cold pools dissipate in less than 30 min.
We see the benefit of CoolDeTA in allowing for a systematic disentanglement of the processes leading up to organized convective cloud and rainfall fields in the systems mentioned above.We have shown that the method can build "cold pool family trees" in RCE and diurnal cycle frameworks, with and without wind shear.These family trees lay the groundwork for a more elaborate analysis of the evolution of convective organization throughout longer simulated periods.Moreover, as CoolDeTA considers dynamic CP signatures, the identified CP boundaries align with the cloud patterns associated with cold pools in satellite images or corresponding simulation output.CoolDeTA thus offers a systematic and objective ground truth labeling for artificial intelligence methods that detect cold pools from simulated cloud fields and that potentially pave the way for future satellite-based CP observations (Hoeller et al., 2024).
As an extension of the current method, it could be interesting to devise a moisture tracking for simulated convective cloud fields.Especially for hysteretic phenomena, such as convective self-aggregation, the long-term memory of mesoscale moisture anomalies could play a critical role in maintaining a persistent atmospheric circulation (Bretherton et al., 2005;Jensen et al., 2021;Muller & Bony, 2015).In RCE, horizontal moisture variations were suggested to give rise to imbalances in thermal radiative emission, which would then drive a persistent circulation and resultant upstream moisture transport-reinforcing the initial moisture imbalance (Emanuel et al., 2014).Recent numerical work extends this notion to diurnal cycle simulations, where MCS emerge and cause abrupt transitions to such a moisture segregated state (Jensen et al., 2021).Importantly, the mesoscale convective system dynamics was found to occur only at sufficiently high horizontal model grid resolution, where CP interactions could be resolved sufficiently.Deciphering the complex moisture and precipitation dynamics requires a tool, such as CoolDeTA, where causal relations can be objectively mapped.

Conclusion
We have devised a CP detection and tracking method, CoolDeTA, which uses a combination of thermodynamic and dynamical variables to track cold pools from their initial forcing by rainfall evaporation until their dissipation stage.The method is shown to function well in a range of contexts, involving RCE and diurnal cycle experiments, as well as simulations with or without wind shear.In comparison to existing buoyancy-focused CP detection methods, CoolDeTA offers the additional benefit of determining the actual dynamical gust front, which is often displaced several kilometers from the thermodynamic boundary.In this way, the method is particularly well suited for the tracking of CP families, where the entire causal cascade of cold pools and their "offspring" can be identified.Such offspring tracking and attribution to parent cold pools would be hampered if only the thermodynamic gust front were available-many new rain cells are triggered at the dynamic gust front.

Figure 1 .
Figure 1.Proposed cold pool (CP) detection and tracking algorithm.(a) Histogram of the indicator field, I seg (i,j,n) (Equation 6), for diu4K at t 471 ."no cold pool (CP)" and "potential CP" grid cells, as determined by the k-means algorithm, are shown as hatched and striped curve fillings.(b) Main algorithm steps: Rain patches (RPs) and CPs are shown as blue and orange areas, respectively.Markers (Ms) are denoted as "x"-symbols; Dashed and dotted lines indicate CP and RP contours, respectively.The steps are repeated for every simulation output time step Δt.(c) Elevation map values (blue curve), I top (AB ↔

Figure 2 .
Figure 2. Comparison of cold pool (CP) detection methods.(a) Time step index 471 of rce0K, showing near-surface virtual temperature perturbations, T′ v , with superimposed dynamical gust front, that is, w > w + 2σ w (red scatter); The superimposed colors represent contours of CP patches ≥A 0 based on different CP detection methods.The enlarged inset shows a magnification for a subregion highlighting a sequence of four time steps.(b) Analogous to (a) but for diu4K.(c) Analogous to (b) but for time step index 485, corresponding to the peak of diurnal rainfall.(d) Analogous to (b) but for time step index 499.(e) Surface temperature, T s , for the two simulation setups visualized in (a)-(d); Solid vertical lines within the plot denote the time steps depicted in (a)-(d).(f) Analogous to (e) but for total domain rainfall per time step, R; The chosen spin-up period of 3 days (432 time steps) is indicated by the vertical dotted line and the horizontal arrow above the panel.

Figure 3 .
Figure 3. Cold pool (CP) expansion for various simulation setups.(a) Median CP area, A cp , versus CP age, t cp for cold pools without any children during their lifetime.Colors represent different simulation setups; Dashed gray lines indicate areas ∝t cp 1.2 ; Only CP ages with ≥50 cold pools are depicted.(b) Analogous to (a) but for cold pools with children during their lifetime.Note the logarithmic vertical axis scaling.

Figure 4 .
Figure 4. Total cold pool (CP) rainfall.(a) Total accumulated rainfall associated with cold pools, R cp , depending on their number of children for different simulation setups.Colored bars represent the interquartile ranges IQR = Q3 Q1 of the simulation setups, with the first quartile Q1 and the third quartile Q3, along with the corresponding median (horizontal dash).Whiskers range from Q1 1.5 × IQR (minimum) to Q3 + 1.5 × IQR (maximum); Outliers w.r.t.this range are not visualized.(b) Analogous to (a) but for the number of parents.Note the logarithmic vertical axis scaling.

Figure 5 .
Figure 5. Children probability distribution.(a) Exceedance probability for a cold pool (CP) to have a certain number of children; Colors indicate different simulation setups.(b) Simple model for CP expansion and spawning of new rain cells.(c) Analogous to (a) but for our simple model.

Figure 6 .
Figure 6.Cascades of tracked cold pools.(a) Labeled cold pools for a subsection of the full simulation domain at time step index 764 of diu4K.The inset plot depicts full domain specific humidity anomalies at the same time step; The black square indicates the subsection of interest.(b) Analogous to (a) but without inset plot and for time step index 772.(c) Analogous to (b) but for time step index 780.(d) Family tree for CP "46" and CP "49."The tree features the IDs of all related cold pools which are present in (a)-(c); Columns represent the time steps indices of CP formation; Filled arrows point from the IDs of CP parents to the IDs of CP offspring; Empty arrows denote additional parents; Red text indicates the IDs of cold pools that are not present in the depicted time steps of (a)-(c), but which contributed to the related cold pools.

Journal of Advances in Modeling Earth Systems
mar ) by determining the nearest integer pixel indices.We collect these marker indices into a field p mar 1d visualizes both cases: On the one hand, RP3 is tracked in time step t n since the corresponding RP overlaps RP3 from the previous time step t n 1 .On the other hand, RP4 is tracked in time step t n since the corresponding RP fully overlaps CP4 from time step t n 1 .If a RP labeled m′ is overlapped sufficiently by a previous RP labeled m 1 and a previous CP labeled m 2 , that is,

with Journal of Advances in Modeling Earth Systems
and if parts of its latest RP p m r from time step t last n are classified as "no CP" in p bin(i, j, n), that is, O seg m,1,n last ,n) < 1.0 × A p m r )HOELLER ET AL.

Table 1
Parameters and Conditions in the Proposed Algorithm Parameters handling rain intensity threshold, minimal potential cold pool (CP) and rain patch (RP) areas, overlap conditions, and CP dissipation.Parameters and conditions below the horizontal line represent additional algorithm rules.
HOELLER ET AL.