Inferring the controlling factors of ice aggregation from targeted cloud seeding experiments

Ice aggregation in clouds plays a crucial role in cloud development and precipitation formation. Despite the significance of ice aggregation, direct in situ quantification of aggregation rates in natural clouds has been challenging due to the difficulty of tracking ice crystals. Here, we present in situ measurements of ice aggregation rates in persistent supercooled stratiform clouds. Using novel glaciogenic seeding experiments (CLOUDLAB), ice crystals are nucleated upwind and subsequently measured downwind after a known advection time in cloud, allowing us to estimate their age. A deep-learning-based detection algorithm (IceDetectNet) counts the individual monomers of aggregates to derive the initial ice crystal number concentration (ICNC${}_{t_{\mathrm{0}}}$). We considered several factors that may influence ice aggregation, including ICNC${}_{t_{\mathrm{0}}}$, temperature, ice crystal size, aspect ratio, and turbulence. Among these, ICNC${}_{t_{\mathrm{0}}}$ was found to be the dominant factor controlling aggregation rates by three independent approaches: causal inference, a physical equation, and machine learning models. We report, however, a subquadratic dependence of the aggregation rate on ICNC${}_{t_{\mathrm{0}}}$ (mean exponent ∼ 0.92; 95 % CI: 0.88–0.97), in contrast to theoretical expectations (quadratic dependence). One possible explanation is that aggregation may also involve smaller ice crystals, but this remains hypothetical. To predict aggregation rates, we evaluated 11 machine learning models and a physically based formulation. CatBoost achieved the best statistical performance, while the physical model proved more robust in sensitivity tests. These findings provide new insights into the microphysical and environmental controls of ice aggregation and establish a robust methodological foundation for studying aggregation processes in natural clouds.

1 Introduction

Ice aggregation is a key microphysical process that influences cloud development and precipitation formation. This process has important implications for weather prediction and climate modeling. During the early stage of ice growth, diffusional processes such as the Wegener–Bergeron–Findeisen mechanism dominate (Korolev, 2007) and vapor deposition in fully glaciated clouds such as cirrus(Gierens et al., 2003). As clouds mature, collisional processes including aggregation and riming become increasingly important (Connolly et al., 2012; Heymsfield, 1986; Sölch and Kärcher, 2011; Hosler et al., 1957). These processes enable ice crystals to grow into larger particles more rapidly than by vapor deposition alone and are fundamental for the formation of snowflakes and graupel, eventually leading to precipitation-sized hydrometeors (Heymsfield, 1986; Sölch and Kärcher, 2011).

Numerous in situ observations have confirmed that ice–ice aggregation occurs in clouds over a wide temperature range, from just below 0 down to −60 °C (e.g., Connolly et al., 2005; Crosier et al., 2011; Field and Heymsfield, 2003). Aggregated ice crystals often constitute a substantial fraction of total ice. For example, irregular ice crystals – which include both aggregated and aged ice crystals – have been reported to account for 84 % of total ice crystals in stratiform clouds (Korolev et al., 2000), 88 % in ground-based snow particle measurements (Zamorsky, 1955), and 94 % in Arctic clouds (Korolev et al., 1999). More specifically, aggregated ice alone has been observed to comprise 38 % of ice crystals in Arctic mixed-phase clouds (Zhang et al., 2024), 45 % in thunderstorm clouds (Jaffeux et al., 2022), and 52 % on average across multiple cloud types during aircraft campaigns (Moss and Johnson, 1994). Despite their apparent ubiquity, quantitative understanding of aggregation rates remains limited.

Several factors have been hypothesized to influence ice aggregation, such as temperature, ice crystal shape, size, and turbulence. Temperature dependence has long been debated, with laboratory studies conflicting trends: Hosler and Hallgren (1960) observed a maximum aggregation efficiency near −15 °C, possibly due to the prevalence of dendritic growth forms that interlock upon collision. In contrast, earlier findings by Hosler et al. (1957) suggested peak aggregation rates around 0 to −5 °C, where quasi-liquid layers on ice surfaces may enhance adhesion (Lamb and Verlinde, 2011). These competing mechanisms – the habit-based, dendritic ice crystal interlocking mechanism at colder temperatures versus the enhanced surface stickiness near freezing, which happens at the same time – highlight the complexity of aggregation processes. Turbulence may further enhance aggregation by introducing small-scale velocity fluctuations that promote collisions (Chellini and Kneifel, 2024; Sheikh et al., 2022). They also suggest that multiple other pathways may operate, depending on the ambient conditions. Ice crystal size and number concentration also play important roles. Larger crystals exhibit greater fall-speed differences and larger geometric cross-sections, while higher ice crystal number concentrations (ICNC) increase collision frequency; both factors enhance the probability of collision and sticking (Hobbs et al., 1974; Field and Heymsfield, 2003; Field et al., 2006; Connolly et al., 2012; Karrer et al., 2021). However, disentangling these effects in natural clouds remains challenging due to observational constraints and the interplay of multiple factors.

Several experimental approaches have been developed to estimate aggregation, yet each carries inherent limitations. Aircraft measurements, often using a “Lagrangian spiral descent” strategy, infer aggregation efficiency from changes in ice crystal size distributions as the aircraft descends at approximately the terminal fall speed of the crystals (Field and Heymsfield, 2003; Field et al., 2006). These measurements rely on 2-D imaging probes and typically assume that decreases in ice crystal number concentration result from aggregation. However, uncertainties arise from experimental artifacts (e.g., ice crystal breaking up on the inlets of the probe) and misattribution of number loss to aggregation alone (McFarquhar et al., 2007; Lawson, 2011). Laboratory studies using ice cloud chambers offer better control over experimental conditions and avoid shattering artifacts, but chamber dimensions are generally too small to accommodate the timescales required for aggregation to occur naturally (Shaw et al., 2020). Moreover, both aircraft and laboratory studies typically rely on indirect estimates – based on changes in ICNC – to infer aggregation efficiency, rather than direct observations (Connolly et al., 2012; Field and Heymsfield, 2003).

These challenges make direct in situ measurements of aggregation rates in natural cloud systems challenging, yet critical for constraining microphysical parameterizations in models. In this study, we use a novel combination of UAV-based glaciogenic cloud seeding experiments (Henneberger et al., 2023) within CLOUDLAB project and an advanced deep-learning detection algorithm (Zhang et al., 2024) to address these gaps. This approach establishes a well-defined initial condition for ice crystal formation and subsequently captures a detailed snapshot of the ice crystal population after their residence in cloud, thereby allowing aggregation to be quantified over a precisely constrained time interval. The deep-learning algorithm (IceDetectNet) robustly identifies individual ice monomers within aggregated ice crystals, providing a monomer-resolved level of detail for quantifying aggregation.

The specific research questions that we address are:

1. Which microphysical and meteorological factors control the rate of ice aggregation?

2. How rapidly does aggregation occur following initial ice formation, and to what extent can the aggregation rate be inferred from the controlling parameters?

To answer these questions, we investigate ice aggregation rates in stratiform clouds using both data-driven and physically derived approaches. Specifically, we (i) identify the microphysical and meteorological controlling controls on aggregation (Sect. 3.1); (ii) disentangle the direct and indirect effects of each factor using causal inference (Sect. 3.2); (iii) evaluate predictive models trained on these factors (Sect. 4.2 and Sect. 4.3); and (iv) test the sensitivity of the predictions to temperature and the initial ice crystal number concentration (ICNC${}_{t_{\mathrm{0}}}$) (Sect. 4.4).

2 The CLOUDLAB campaign and observational data

The CLOUDLAB experiments were conducted in persistent wintertime stratus clouds over the Swiss Plateau, with measurements centered at a site near Eriswil (47°04^′14^′′ N, 7°52^′22^′′ E; 920 m a.s.l.). These clouds were typically supercooled, liquid-dominated, and quasi-stationary, with bases below 1000 m a.g.l. and thicknesses of several hundred meters (Scherrer and Appenzeller, 2014). We performed glaciogenic seeding upwind of the measurement site using uncrewed aerial vehicles (UAVs) and sampled the resulting microphysical changes with a suite of in situ and remote sensing instruments. Methods specific to this study are detailed below, including the seeding operations and instrumentation (Sect. 2.1) and observed ice properties (Sect. 2.4).

2.1 Seeding Operations and Instrumentation

Glaciogenic seeding was performed upwind of the measurement site using a customized uncrewed aerial vehicle (UAV; Meteodrone MM-670, Meteomatics AG, Switzerland) equipped with burn-in-place flares (Zeus MK2, Cloud Seeding Technologies, Germany). Each flare contained approximately 200 g of seeding material, including about 20 g of silver iodide and other ice-active compounds effective at temperatures below −5 °C (Chen et al., 2024; Miller et al., 2024). The UAV was operated by flying multiple crosswind legs (200–400 m) while releasing seeding particles for 5–6 min at distances of 2–3 km upwind of the measurement site at seeding time t₀ (In-cloud seeding and ice nucleation processes in Fig. 1). Further details on the seeding operations and ice nucleation mechanisms are provided in Miller et al. (2024, 2025).

Figure 1 Scheme for quantifying in situ ice aggregation after cloud seeding. Seeding particles are released by a UAV, initiating ice crystal formation at time t₀. The ice crystals grow and aggregate during the advection time t₁−t₀, and at t₁ an image of them is captured by a holographic imager on a tethered balloon system. The ice habit and the number of monomers per aggregate are quantified based on the IceDetectNet (Zhang et al., 2024). For example, the right panel shows an aggregate consisting of two monomers classified as a frozen droplets and a column-plate. The aggregation rate (R_agg) is defined as the total number of aggregation events, $\sum (n_i-\mathrm{1})$, divided by the advection time t and the sampled volume $V_{t_{\mathrm{1}}}$ at t₁.

The plume of ice crystals generated by seeding arrived at the measurement site after 5–10 min depending on the wind speed. Because the plume was typically several hundred meters wide and persisted over the site for several minutes, it could be sampled continuously as it passed, yielding a time series of microphysical observations. The advection time t – the time between ice nucleation and observation – was calculated by comparing the times of seeding particles release and first ice detection, and using the observed wind speed across the plume to estimate advection time. To avoid underestimating travel time, the maximum wind speed among the measured instruments was used in the calculation (see Fuchs et al., 2025 for more details regarding the method). Throughout this growth period, the temperature and wind speed were assumed to remain constant.

The microphysical properties of the plume were measured using a suite of in situ and remote sensing instruments at measurement time t₁ (see Fig. 1, Ice detection). Remote sensing measurements included ground-based Ka- and W-band Doppler cloud radars. These radars were used to observe the general cloud structure (e.g., cloud top height) before and after seeding.They were also used to retrieve the turbulence intensity, which is expressed as the eddy dissipation rate (EDR). EDR was retrieved from the Ka- and W-band radars following the spectral-width method (Appendix E in Wu et al., 2025). Primary estimates were obtained from the Mira35 MBR7 Ka-band radar. In addition, the RPG94 W-band radar was used to supplement data gaps. The retrieval was temporally averaged over 30 s, and the resulting EDR fields have a spatial resolution of approximately 30 m.

In situ observations were provided by a tethered balloon system (TBS) carrying the HOLographic Imager for Microscopic Observations (HOLIMO; Ramelli et al., 2020). The TBS consisted of a 200 m³ helium-filled kytoon capable of reaching altitudes up to 1 km above ground. Suspended 30 m below the kytoon, the instrument platform carrying HOLIMO, a digital in-line holography system with resolving sizes of 6 µm for cloud droplets and 25 µm for ice crystals. HOLIMO operated continuously before, during, and after each seeding experiment, recording cloud droplets throughout and ice crystals over the duration of the plume encounter. Holograms were reconstructed using HoloSuite (Fugal and Shaw, 2009) at 1 s resolution, and particles were classified as artifacts, cloud droplets, or ice crystals using a convolutional neural network (Touloupas et al., 2020). Manual verification was performed for all cloud droplets and ice crystals larger than 35 µm to ensure classification accuracy; no ice crystals smaller than this threshold were observed.

Classified ice crystals were then further analyzed using a fine-tuned version of IceDetectNet (Zhang et al., 2024), named IceDetectNet-CLOUDLAB (see Sect. 2.3 for details). IceDetectNet-CLOUDLAB is a rotated object detection model trained to assign each monomer a shape label comprising its basic habit and microphysical process, and to estimate the number of monomers per ice crystal. The classification scheme includes three basic ice habits – column, plate, and irregular – and one microphysical process, riming. If an ice crystal contains more than one monomer, it is additionally classified as aggregated. This results in two independent microphysical processes (riming and aggregation) and four possible states for each ice crystal: pristine, rimed, aggregated, and rimed + aggregated. Combining the three habits with these four states yields a total of 12 distinct ice crystal classes. The uncertainty in cloud droplet number concentration is approximately ±5 %, while that for ice crystal number concentration ranges from 5 %–10 % for crystals larger than 100 µm and about 15 % for smaller ice crystals. Uncertainty quantification for habit/process classification and monomer counts is detailed in Appendix A.

2.2 Estimating Aggregation Rates from Experimental Data

Ice aggregation was quantified by estimating the number of monomer-level collisions that formed each observed aggregate. Under the definition that all ice crystals are monomers at the initial time t₀ and that aggregation proceeds via binary collisions, an observed aggregate containing n_i monomers has undergone (n_i−1) aggregation events to reach its observed state at time t₁. The aggregation rate R_agg (s⁻¹ L⁻¹) was thus defined as the mean number of aggregation events per unit time per unit cloud volume, averaged over the advection time t:

$$\begin{array}{}\text{(1)}& R_{\text{agg}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(n_i-\mathrm{1})}{t\cdot V_{t_{\mathrm{1}}}}}\end{array}$$

where N is the total number of detected aggregates, n_i is the monomer count of the ith aggregate, t is the advection time estimated based on the method in (Fuchs et al., 2025), and $V_{t_{\mathrm{1}}}$ is the sampled cloud volume at time t₁. The approach relies on the following assumptions:

1. Steady background conditions. Environmental and microphysical properties of the background cloud (e.g., temperature, wind speed) are assumed to remain constant throughout each experiment.

2. Instantaneous nucleation. Ice nucleation is assumed to occur immediately upon release of seeding particles, providing a well-defined initial time to ice formation. Our observations have revealed that ice formation is likely initiated very quickly after the release of the highly hygroscopic and ice-active seeding material (Miller et al., 2025).

3. Conservation of monomer number. The total number of ice crystal monomers is assumed to be conserved between the seeding and observation locations. This implies that (i) sedimentation losses within the plume are approximately balanced by ice crystals falling from above. This assumption is supported by radar observations, which do not show a systematic downward displacement of the enhanced reflectivity associated with the seeded plume during the observation period (see Fig. 2 in Fuchs et al., 2025). (ii) no secondary ice production (SIP) occurs. The latter is supported by observations: no large droplets or graupel were detected, and no evidence of frozen droplet breakup was observed, consistent with conditions unfavorable for SIP.

4. Complete detection of aggregate monomers. All monomers within an aggregate are assumed to be accurately identified by IceDetectNet. However, some early aggregates may become unresolved after diffusional growth and be classified as irregular, but this fraction is very small (see Fig. 2h). Other potential detection errors and classification errors are discussed in Appendix A and are considered negligible.

2.3 Training IceDetectNet-CLOUDLAB for Aggregated Ice Monomer Identification

To retrieve the number of monomers in each aggregate, we fine-tuned the original IceDetectNet model (Zhang et al., 2024) to create IceDetectNet-CLOUDLAB, specifically adapted for holographic images collected during CLOUDLAB seeding experiments. The original model was trained to classify ten ice crystal habits; for this application, the output layer was reconfigured to distinguish four categories: column, plate, column-rimed, and plate-rimed. The final classification layer was randomly initialized, while all other parameters were taken from the IceDetectNet. A total of 2380 manually labeled images were used, randomly selected to ensure an approximate balance across all seeding experiments. Of these, 2134 images were used for training and 246 for testing. The model architecture followed IceDetectNet-CLOUDLAB (Zhang et al., 2024), based on S2ANet with a ResNet-50 backbone. To improve training stability on the smaller CLOUDLAB dataset, which contains fewer ice crystals and fewer classes than previous applications, the initial learning rate was reduced to 0.0001. The number of training epochs was extended to 200, with learning rate decay scheduled earlier at epochs 32 and 48 to encourage earlier convergence and to prevent overfitting. A linear warmup phase of 1000 iterations was applied to avoid early gradient instability. During inference, up to six predictions per image were retained after non-maximum suppression (IoU threshold = 0.5).

2.4 Statistical Characterization of Ice Crystal Properties

We analyzed 21 seeding experiments (Table F1) conducted at temperatures between −4.7 and −7.8 °C. The ice crystal habit distribution differed systematically with temperature: warmer experiments (T> −7 °C) contained exclusively columnar crystals, whereas the colder experiments (T≤ −7 °C) contained both plates and short columns (Fig. E1). Plates and columns have been shown to exhibit different diffusional growth rates along their major axes (Fuchs et al., 2025), which can influence subsequent aggregation. Based on this habit difference, the experiments were classified into “warmer” (T> −7 °C; 14 experiments, 2508 data points) and “colder” (T≤ −7 °C; 7 experiments, 797 data points) regimes (Fig. 2a). Advection times were generally shorter in the colder experiments, ranging from 294–532 s, compared to 371–629 s in the warmer experiments (Fig. 2b).

ICNCt₀, the initial ice crystal number concentration, was estimated by assuming that each aggregated monomer corresponds to one ice crystal present at t₀. This allows estimation of the ICNC${}_{t_{\mathrm{0}}}$ at the initial state t₀. ICNC${}_{t_{\mathrm{0}}}$ was higher in the warmer regime than in the colder, with mean ± standard deviation of approximately 815 L⁻¹ ± 394 L⁻¹ and 566 L⁻¹ ± 291 L⁻¹, respectively (Fig. 2c). While these concentrations are higher than those typically found in unseeded clouds (Heymsfield and Willis, 2014), the aggregation process is governed by the same underlying collisional physics.Aggregates were more frequent in the warmer regime, with fractions reaching up to 43 % of total observed ice crystals, while remaining around 15 % in the colder regime (Fig. 2d). Ice crystals in the warmer regime reach mean major size (defined as the length of the crystal’s major axis) per experiment of up to  312 µm (minimum  101 µm), compared to a range of  66–129 µm in the colder regime; the warmer regime also exhibited a broader, long-tailed size distribution (Fig. 2e). Aspect ratio (defined as the ratio of the major to minor axis lengths) distributions differed as well: the colder regime was narrowly peaked around 1.4–1.6, while the warmer regime showed a broader and flatter distribution around 1.73–4.13 (Fig. 2f). The riming ratio (RR) is defined as the fraction of rimed crystals. Rimed crystals were abundant in both regimes, with fractions exceeding 38 % in all experiments and reaching up to 94 % in one single experiment, indicating that riming was nearly ubiquitous (Fig. 2g). Finally, irregular ice crystals, which have no identifiable habit, were uncommon overall (2 %–5 %; Fig. 2h). Some early aggregates may have grown sufficiently by diffusion to obscure their monomer structure, thus being classified as irregular. These crystals were slightly more frequent in the colder regime, suggesting that any resulting underestimation of aggregation is likely negligible. EDR values measured during the seeding flights were modest (on the order of 10⁻⁴–${\mathrm{10}}^{-\mathrm{3}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}}$, Fig. B1), comparable to turbulence levels reported for other boundary-layer cloud environments (Chu et al., 2025). This indicates that aggregation developed under a weak-turbulence conditions characteristic of stratiform mixed-phase clouds.

Figure 2 Distributions of ice crystal properties during seeded periods, grouped by temperature: colder (T≤ −7 °C, blue) and warmer (T> −7 °C, red). (a) Temperature; (b) Advection time; (c) ICNC${}_{t_{\mathrm{0}}}$; (d) Aggregation ratio; (e) Major size of ice crystals; (f) Aspect ratio; (g) Riming ratio; (h) Irregular ratio. For (a) and (b), the solid lines show the per-experiment distributions. For (c) through (h), the solid lines represent one-second resolution distributions within each regime, and the scatter points indicate the means of the experiments (The scatter positions along the y axis are offset solely to reduce overlap. Only the values along the x axis carry meaning).

3 Microphysical and Environmental Controls on Ice Aggregation

We examined the microphysical and environmental controls on ice aggregation through a structured sequence of analyses. We first examined how microphysical and environmental factors influence the observed aggregation rates (Sect. 3.1). We then applied a causal graph to disentangle the direct and indirect effects of these factors on aggregation (Sect. 3.2).

3.1 Correlations Between Aggregation Rate and Microphysical and Environmental Factors

To better understand the mechanisms underlying aggregation, we first examine the Pearson correlations between the aggregation rate and four factors: ICNC${}_{t_{\mathrm{0}}}$, advection time, temperature, and EDR. No clear correlation was observed between aggregation rate and EDR, which characterizes the turbulence intensity within the cloud (see Appendix  B for the EDR correlation analysis). This suggests that, under the present experimental conditions, turbulence at the resolved scales (30 m × 30 m) did not significantly influence aggregation. The correlations with ICNC${}_{t_{\mathrm{0}}}$, advection time, and temperature are all presented in the following.

3.1.1 Strong positive correlation between ICNC${}_{t_{\mathrm{0}}}$ and aggregation rate

ICNC${}_{t_{\mathrm{0}}}$ exhibited a strong positive correlation with the aggregation rate in all experiments (Fig. D1). This is consistent with theoretical expectations, as higher crystal concentrations increase the probability of collisions leading to aggregation (Hobbs et al., 1974).

To illustrate this relationship more clearly, we present two experiments for reference: a warmer case (SM054, T= −5 °C) and a colder case (SM069, T= −7.6 °C) (Fig. 3a, b). Under the condition that the ice crystals are homogeneously distributed, the instantaneous aggregation rate (Seifert and Beheng, 2006, Eq. 62), which is the number of aggregation events per unit time per unit volume of cloud, can be expressed as:

$$R_{\text{agg}}=\mathrm{1}/\mathrm{2}K(D_{\mathrm{1}},D_{\mathrm{2}},T)N_{\mathrm{1}}{N}_{\mathrm{2}},$$

where N₁ and N₂ are the number concentrations of two ice crystal populations of sizes D₁ and D₂, T is temperature, and $K(D_{\mathrm{1}},D_{\mathrm{2}},T)$ is the collision kernel, which depends on ice crystal size, shape, fall speed, and temperature (Connolly et al., 2012). In the simplest case, if the shape of the size distribution is approximately preserved and only its overall magnitude changes with the total ice number concentration, then the concentrations of the relevant size classes, N₁ and N₂, both scale linearly with the initial total ice concentration, ICNC${}_{t_{\mathrm{0}}}$. Under this approximation, both the instantaneous and time-averaged aggregation rates are expected to scale as $R_{\text{agg}}=\mathrm{1}/\mathrm{2}K(D_{\mathrm{1}},D_{\mathrm{2}},T)({\text{ICNC}}_{t_{\mathrm{0}}}{)}^{\mathrm{2}}$.

Figure 3 Correlation between aggregation rate and ICNC${}_{t_{\mathrm{0}}}$ in two temperature regimes. (a) warmer case (SM054, T= −5.0 °C); (b) colder case (SM069, T= −7.6 °C). One-second data points are shown as red circles (warmer) and blue squares (colder). Solid lines with shaded areas show free power-law fits $R_{\text{agg}}=a\cdot {\text{ICNC}}_{t_{\mathrm{0}}}^n$ with 95 % confidence intervals (red for warmer, blue for colder); dashed lines with grey shading show quadratic fits (n=2). The prefactor a was fitted independently for each experiment and reflects the effective collision kernel.

Motivated by this scaling, we first evaluated the data against this fixed quadratic relationship (R_agg $=a\cdot$ ICNC${}_{t_{\mathrm{0}}}^{\mathrm{2}}$), where the prefactor a effectively represents $\mathrm{1}/\mathrm{2}K(D_{\mathrm{1}},D_{\mathrm{2}},T)$. To account for differences between experiments, a was fitted independently in each case. While this model gets correlation coefficients of r=0 and r=0 for the warmer and colder cases, respectively, it systematically underestimated aggregation at low ICNC${}_{t_{\mathrm{0}}}$ (Fig. 3a, b, dashed grey lines).

A better fit was achieved using a free power-law relationship (R_agg $=a\cdot$ ICNC${}_{t_{\mathrm{0}}}^n$), where both a and n were estimated separately for each experiment. This empirical model closely matched the observations across the full range of concentrations (Fig. 3a, b, solid lines). In the warmer-case experiment (SM054), the aggregation rate scaled with an exponent of 1, improving the correlation to r=0 (Fig. 3a). In the colder-case experiment (SM069), the best-fit exponent was 0, with a corresponding correlation of r=0 (Fig. 3b). This pattern was consistent across all 21 experiments, where the free-fit power-law systematically performed better than the fixed-quadratic form (see Fig. D1). The average exponent across all experiments was 0 ± 0, smaller than the quadratic dependence expected from the collision kernel formulation (Melzak, 1957; Connolly et al., 2012; Seifert and Beheng, 2006). Bulk microphysics schemes such as Lin et al. (1983) and Morrison and Milbrandt (2015) have even adopted a linear dependence on ICNC, which is closer to our observations, often with an additional size threshold so that aggregation only occurs once ice crystals exceed a certain size. However, our observations hint that aggregation may also occur among smaller ice crystals. This could result in effective dependence with an exponent smaller than one.

3.1.2 No significant correlation between aggregation rate and advection time

Advection time in our experiments ranged from 294 to 650 s (Table F1, Fig. 4) and showed no significant correlation with aggregation rate (r= −0.09, p=0). This parameter is specific to our experimental design, as the ice crystals were artificially generated by cloud seeding, enabling an explicit measurement of the total growth period. In natural-cloud studies, advection time is rarely quantified, and direct comparisons are therefore not straightforward. Theoretically, a longer advection time would be expected to yield a smaller time-averaged aggregation rate, because ICNC decreases over time. We find only a very weak tendency for R_agg to decrease with advection time (Fig. 4, R²=0), suggesting that our results remain broadly comparable to previous observational and laboratory studies despite this experimental uniqueness (Connolly et al., 2012; Field and Heymsfield, 2003; Field et al., 2006). This weak dependence between aggregation rate and advection time may reflect the fact that ice crystals must first grow to sufficiently large sizes before collisional aggregation becomes efficient. During the early growth phase, ICNC decreases rapidly due to plume dilution (as seeded particles are initially concentrated and then disperse; see Fig. 2 in Ramelli et al., 2024), but this shows little impact on the aggregation rate, making it appear largely insensitive to advection time.

Figure 4 Correlation between aggregation rate and advection time. Orange markers with error bars show the mean ± standard deviation for each experiment, and the solid line represents the linear fit (r=0, p=0) across all 21 experiments.

3.1.3 Temperature dependence of aggregation rate associated with ice crystal size and aspect ratio

To investigate the role of temperature in shaping aggregation, motivated by the higher aggregation rate observed in the warmer case (SM054) compared to the colder case (SM069), we first examined how ice crystal size and geometry differ between warmer and colder experiments. Using the same representative cases shown in Fig. 3, we compared the warmer case SM054 (−5 °C) with the colder case SM069 (−7.6 °C). Consistent with known temperature-dependent ice habits, the warmer group (above −7 °C) consisted almost exclusively of columnar ice crystals (Fig. E1a), whereas the colder group (below −7 °C) included both columnar and plate-like crystals (Fig. E1b). Correspondingly, the warmer case exhibited consistently broader distributions of both major axis length and aspect ratio (Fig. 6). This behavior was not unique to these two cases: across the full dataset, warmer experiments showed broader distributions of major size (Fig. 2e) and aspect ratio (Fig. 2f) than colder experiments. Increased variability in crystal size and shape broadens the range of fall velocities among ice crystals (Heymsfield, 1972; Mitchell, 1996), which is expected to enhance collision frequency and thereby promote aggregation.

Consistent with this interpretation, aggregation rate shows a very weak positive association with temperature across all experiments (Pearson r=0, p=0; Fig. 5). While this relationship does not reach conventional significance levels and may therefore be difficult to distinguish from noise over the limited temperature range sampled, it is qualitatively consistent with the observed structural differences in the ice population and with previous laboratory studies suggesting enhanced aggregation under warmer conditions (Hosler et al., 1957).

Figure 5 Correlation between aggregation rate and temperature. Same as Fig. 3, but for the Temperature.

Figure 6 Distributions of ice crystal major axis length and aspect ratio in warmer and colder cases. (a) Major axis length and (b) aspect ratio for SM054 (warmer, red) and SM069 (colder, blue).

3.2 Causal Pathways Among Aggregation Drivers

Previous analyses showed that ICNC${}_{t_{\mathrm{0}}}$, temperature, major axis length, and aspect ratio each correlate with aggregation rate to varying degrees. However, these factors are not independent, and their relative contributions and interactions remain unclear. To disentangle direct and indirect effects among these variables, we constructed a causal graph in the form of a directed acyclic graph (DAG; Fig. 7), incorporating these 4 variables as well as riming ratio (RR), which is the percentage of rimed ice crystals to all the observed ice crystals. Because riming influences ice crystal geometry – particularly major size and aspect ratio – and it may also influence surface stickiness, especially shortly after the riming event, it can, in principle, affect aggregation indirectly.

DAGs represent hypothesized cause–effect relationships, with nodes denoting variables and arrows indicating the direction of influence (Pearl, 2009; Peters et al., 2016). Unlike correlation-based analyses, DAGs explicitly separate direct effects from mediated pathways, enabling quantitative decomposition of total effects. For example, temperature may influence aggregation rate both directly and indirectly by altering ice crystal size and aspect ratio.

The graph structure was specified by combining prior physical knowledge with statistical dependencies inferred from the data. Specifically, candidate causal links were first proposed based on established microphysical understanding of ice aggregation (e.g., ICNC, size, and habit influencing collision and sticking probability), while arrows implying physically implausible causality (e.g., aggregation rate influencing temperature) were excluded a priori. Statistical associations were then used to assess whether the proposed links were supported by the data. The final DAG therefore reflects a physically constrained causal hypothesis rather than a data-driven structure learned automatically. Path coefficients were estimated using structural equation modeling (SEM), which fits a system of regression equations to quantify the strength of each link. All variables were standardized prior to fitting, so that each coefficient represents the change in aggregation rate (in standard deviations) associated with a one-standard-deviation increase in the predictor, reflecting their relative importance.

The fitted DAG (Fig. 7) shows that ICNC${}_{t_{\mathrm{0}}}$ (+0.81), temperature (+0.24), and major size (+0.08) have positive direct effects on aggregation rate, with ICNC${}_{t_{\mathrm{0}}}$ clearly dominating. Aspect ratio exhibits a weak negative direct effect (−0.16), likely because stronger aggregation tends to increase the number of monomers per particle, producing more compact and rounded aggregates and thus reducing the aspect ratio. Beyond its direct effect, temperature also strongly influences aspect ratio (+0.61) and major size (+0.40), which in turn affect aggregation. Since aspect ratio and major size have opposing effects (negative and positive, respectively), these indirect contributions partially cancel, and temperature maintains an overall positive influence.

Figure 7 Causal graph of factors influencing aggregation rate. The graph depicts the inferred direct and indirect effects on aggregation rate (R_agg) among five variables: ICNC${}_{t_{\mathrm{0}}}$, temperature (T), major axis length (Majsiz), aspect ratio (AR), and riming ratio (RR). Nodes represent variables, and directed arrows indicate causal relationships. Arrow thickness, color, and the numerical labels on the arrows denote the magnitude and direction of the standardized effects (change in aggregation rate in standard deviations per standard deviation change in the predictor).

RR has a negligible direct effect on aggregation (−0.01). RR indirectly increases the major axis length (+0.31), which tends to promote aggregation; however, the positive contribution of major size to the aggregation rate (+0.08) is itself very small. Moreover, most rimed ice crystals in our measurements exhibit only light riming (Fig. E1), limiting the potential for riming-induced enhancements in fall speed or sticking efficiency. Under these specific experimental conditions, the overall influence of RR on aggregation appears minimal, suggesting that any coupling between riming and aggregation is weak or not detectable within these experiments.

4 Aggregation Rate Predictive Models

Building on the preceding analysis, we next evaluated whether aggregation rates can be quantitatively predicted from the identified microphysical and meteorological factors. We first describe the training procedures and evaluation metrics (Sect. 4.1). Then, we implemented and compared eleven machine learning models (Sect. 4.2) and one physically based equation (Sect. 4.3), using temperature, major axis length, aspect ratio, and ICNC${}_{t_{\mathrm{0}}}$ as input features and the aggregation rate as the target variable. RR was excluded from the predictive modeling as its direct effect on aggregation was found to be negligible in our causal inference analysis (Sect.3.2). Finally, to evaluate the robustness of these models under varying atmospheric and microphysical conditions, we performed a sensitivity analysis with respect to temperature and ICNC${}_{t_{\mathrm{0}}}$ (Sect. 4.4).

4.1 Predictive Model Training and Evaluation

We trained the models using temperature, major axis length, aspect ratio, and ICNC${}_{t_{\mathrm{0}}}$ as predictors of the aggregation rate. All variables except temperature were available at 1 s resolution; temperature was constant for each experiment. We fitted the analytical model (Eq. 2) by minimizing the mean absolute error (L1 norm) between predicted and observed aggregation rates, with L₂ regularization applied to all parameters (penalty coefficient $={\mathrm{10}}^{-\mathrm{4}}$). All coefficients were initialized to one and optimized using gradient descent with a learning rate of 10⁻² over 1500 iterations. All machine learning models were trained with default hyperparameters, except for a fixed random seed (42) to ensure reproducibility. No early stopping, feature engineering, or categorical coding was applied, as all inputs were continuous. Model performance was assessed using five-fold cross-validation with group-based splits defined by seeding experiment identifiers. This ensured that data from a given experiment were restricted to either training or test sets. For each model and fold, we computed the root mean square error (RMSE), which quantifies the typical magnitude of prediction error, and the coefficient of determination (R²), which reflects the proportion of variance explained. We reported the mean and standard deviation of both metrics across the five folds. After evaluation, each model was retrained on the full dataset and archived for reproducibility. Model interpretability was assessed using SHapley Additive exPlanations (SHAP; Lundberg and Lee, 2017). SHAP values quantify how much each feature shifts a prediction away from the dataset’s average prediction, with positive values indicating an increase and negative values a decrease. Conceptually, each feature is treated as a “player” in a cooperative game, and its SHAP value represents the average change in the prediction when the feature is added to all possible subsets of other features. This provides a consistent measure of each feature’s contribution to the prediction, accounting for interactions with other features. Global feature importance was obtained by averaging the absolute SHAP values across all samples, and the results were visualized using summary plots.

4.2 Machine Learning Models for Aggregation Rate Prediction

4.2.1 Overview of Machine Learning Algorithms

We implemented eleven supervised regression algorithms to model the relationship between environmental and microphysical predictors and aggregation rate. These included linear methods (Linear, Ridge, and Bayesian Ridge Regression), instance-based learning (K-Nearest Neighbors), decision tree–based models (Decision Tree, Extremely Randomized Trees), and ensemble boosting techniques (Adaptive Boosting, Gradient Boosting, Light Gradient Boosting Machine, Extreme Gradient Boosting, and gradient boosting with categorical features support (CatBoost)). Linear models assume independent, additive effects and are limited in representing nonlinear interactions. Tree-based models capture nonlinear and interaction effects through hierarchical partitioning of the feature space, while boosting methods iteratively refine predictions by combining multiple weak learners. Detailed descriptions of all models are provided in Appendix C. Among these, CatBoost is particularly well-suited for the current problem because of its ability to handle nonlinear interactions and its robustness to overfitting. CatBoost is a gradient boosting method that uses symmetric (oblivious) decision trees, where all nodes at the same depth split on the same feature and threshold. This symmetry simplifies optimization and improves generalization. It also introduces “ordered boosting,” which builds trees in a way that avoids using future information during training, thereby reducing overfitting. Additionally, CatBoost efficiently handles categorical features and typically requires minimal hyperparameter tuning.

4.2.2 Predictive Performance and Feature Contributions

Among all machine learning models, CatBoost and Gradient Boosting achieved the highest performance, with the lowest RMSE ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{1}}$ s⁻¹ L⁻¹) and high R² (0) for CatBoost, and a slightly higher RMSE but marginally better R² (0) for Gradient Boosting. Their comparable skill likely reflects their shared ensemble structure and residual learning strategy. Extra Trees and Light Gradient Boosting Machine followed closely in performance (Table 1). For subsequent analyses, CatBoost was selected as the representative data-driven model owing to its competitive accuracy and robust generalization. The physically derived equation also performed comparably to the tree-based models, achieving an RMSE of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{1}}$ s⁻¹ L⁻¹ and R² of 0, suggesting that it captures the dominant aggregation dependencies despite its reduced complexity and lack of data-driven tuning. In contrast, linear models (Linear Regression, Ridge, Bayesian Ridge) performed poorly, with R² values around 0 and RMSE around $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{1}}$ s⁻¹ L⁻¹, highlighting the importance of nonlinear interactions and feature dependencies that linear approaches cannot represent. In each experiment (Fig. 8), the agreement between observed aggregation rates and CatBoost predictions is evident, with residuals having means near zero and relatively small standard deviations. These results demonstrate that the model effectively captures the relevant physical mechanisms.

Table 1 Model performance comparison. Mean RMSE (10⁻¹ s⁻¹ L⁻¹) and R² values along with their standard deviations, sorted by RMSE in ascending order.

Figure 8 Observed, predicted, and residual aggregation rates per experiment. For each experiment, observed aggregation rates (dark blue circles), CatBoost-predicted aggregation rates (yellow triangles), and residuals (pink squares, observed minus predicted) are shown side by side, with observed values on the left, predicted in the middle and residuals on the right. Markers indicate mean values, and error bars represent the standard deviation within each experiment. Rates are expressed in L⁻¹ s⁻¹.

To interpret the predictions, we computed SHAP, which quantifies the relative contribution of each feature factor. ICNC${}_{t_{\mathrm{0}}}$ was the dominant predictor, explaining 69.9 % of the variance (Fig. 9), consistent with the earlier correlation analysis (Figs. 3, D1) and the causal graph (Fig. 7), where ICNC${}_{t_{\mathrm{0}}}$ also showed the strongest standardized effect. However, a simple linear regression model using ICNC${}_{t_{\mathrm{0}}}$ alone performed poorly ($R^{\mathrm{2}}=\mathrm{0}\pm \mathrm{0}$, RMSE = $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}\pm \mathrm{0}\times {\mathrm{10}}^{-\mathrm{1}}$ s⁻¹ L⁻¹; Table 1), underscoring the importance of nonlinear interactions among multiple variables.

Figure 9 SHAP analysis of feature contributions to CatBoost model predictions. Each dot represents one prediction at 1 s resolution from the test data. SHAP values (horizontal axis) indicate the impact of each input feature on the predicted aggregation rate. Features (i.e. ICNC${}_{t_{\mathrm{0}}}$, Major size, Aspect ratio, and Temperature) are ranked by their overall importance (right), with percentages representing the mean absolute SHAP value normalized across all features. Dot colors reflect the corresponding feature values, with red denoting higher and blue denoting lower values of the feature.

Major axis length (14.9 %) also contributed meaningfully to the predictions. Larger major sizes were associated with higher aggregation rates, consistent with the expectation that larger crystals, with greater cross-sectional area and fall-speed variability, enhance collision likelihood (Heymsfield and Miloshevich, 2003) – in agreement with the positive direct effect identified in the causal graph (+0.08, Fig. 7). Aspect ratio (8.3 %) had a weaker and more complex effect: SHAP values clustered near zero, with a slight tendency for lower aspect ratios to favor aggregation and higher aspect ratios to suppress it, consistent with its negative effect in the causal graph (−0.16).

Temperature accounted for 6.9 % of the model output. Higher temperatures generally promoted aggregation, consistent with laboratory findings (Hosler et al., 1957) and the causal graph (+0.24, Fig. 7). Also, unlike the other variables, temperature was measured as a single value per experiment rather than at a 1 s resolution, which may have limited its explanatory power in the model. Nevertheless, its limited contribution suggests that its influence is largely indirect and mediated by changes in crystal habit and size, as reflected by aspect ratio and major axis length. Therefore, temperature likely acts as a secondary driver embedded within other structural parameters.

4.3 Physical Model

4.3.1 Model Formulation

We constructed an empirical formulation for the ice aggregation rate based on established microphysical principles. The basic structure follows the collision kernel framework (Melzak, 1957), in which the aggregation rate is proportional to the product of the number concentrations of the interacting particles. Under the simplest assumption of equal-sized populations, this scaling becomes quadratic in ICNC${}_{t_{\mathrm{0}}}$. However, several bulk microphysics schemes adopt a linear dependence. For example, in Lin et al. (1983), the aggregation rate is parameterized as being proportional to the ice mixing ratio (Lin et al., 1983, Eq. 21). Because the mixing ratio is given by the product of number concentration and mean ice crystal mass, this formulation effectively yields a linear dependence on number concentration when the ice-crystal-mass–number relationship (i.e., the assumed PSD shape) is held fixed within the scheme. Similar linear forms also appear in two-moment parameterizations such as Morrison and Milbrandt (2015). To allow flexibility between these theoretical and parameterization-based scalings, the exponent on ICNC${}_{t_{\mathrm{0}}}$ was treated as a free parameter. The temperature dependence is represented as an exponential term, consistent with parameterizations of collision efficiency in two-moment schemes (e.g. Seifert and Beheng, 2006), where temperature modulates the quasi-liquid layer and hence the sticking probability upon collision. The dependence on major axis length is formulated as a power law to account for its role in determining both geometric cross-section and differential fall speed, which are classical components of aggregation kernels. Thus, the proposed equation is:

$$\begin{array}{}\text{(2)}& R_{\text{agg}}=\mathit{\alpha }\cdot {\mathrm{ICNC}}_{t_{\mathrm{0}}}^{{\mathit{\beta }}_{\mathrm{0}}}\cdot {\exp }^{{\mathit{\beta }}_{\mathrm{1}}T}\cdot {\text{MajSiz}}^{{\mathit{\beta }}_{\mathrm{2}}}\end{array}$$

where: R_agg is the observed aggregation rate (s⁻¹ L⁻¹), α is a dimensioned scaling constant, β₀, β₁, and β₂ are empirically fitted exponents, ICNC${}_{t_{\mathrm{0}}}$ is the initial ice crystal number concentration (L⁻¹), T is the ambient temperature (°C), MajSiz is the average major axis length of ice crystals (m). The term ICNC${}_{t_{\mathrm{0}}}^{{\mathit{\beta }}_{\mathrm{0}}}$ accounts for the increased probability of collisions as a function of ice crystal concentration. Collection theory suggests a near-quadratic dependence under the assumptions of random motion and homogeneous mixing (Connolly et al., 2012). We retain a flexible exponent β₀ to account for inhomogeneities and dispersion effects in real clouds. The temperature-dependent exponential term, $e^{{\mathit{\beta }}_{\mathrm{1}}T}$, reflects the nonlinear role of temperature in aggregation-relevant processes. Temperature modulates diffusional growth, growing dimension and the thickness of the quasi-liquid layer, which shapes ice size and habit. Many of these processes exhibit exponential or Arrhenius-like behavior with temperature. The exponential formulation is consistent with aggregation parameterizations in operational two-moment microphysics schemes (e.g. Seifert and Beheng, 2006; Lin et al., 1983). The power-law dependence on ice crystal size, MajSiz${}^{{\mathit{\beta }}_{\mathrm{2}}}$, represents the combined influence of fall speed variability and geometric collision cross section. Larger crystals sediment faster and offer larger interaction surfaces, increasing the likelihood of collision and adhesion. This formulation follows classical approaches to the collection kernel under differential sedimentation. Aspect ratio was excluded from the final model, as it is not explicitly represented in two-moment bulk schemes (Lohmann and Roeckner, 1996; Seifert and Beheng, 2006), and its fitted coefficient in our analysis was negligible compared to the other factors.

The final fitted form is:

$$\begin{array}{}\text{(3)}& R_{\text{agg}}=\mathrm{0}\cdot {\mathrm{ICNC}}_{t_{\mathrm{0}}}^{\mathrm{0}}\cdot \exp (\mathrm{0}\cdot T)\cdot {\text{MajSiz}}^{\mathrm{0}},\end{array}$$

When optimized on the same dataset, this physical model achieved an RMSE of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ L⁻¹ and an R² of 0.78 (Table 1). Although lower than the CatBoost model (R²=0), its performance remains strong, suggesting that the dominant dependencies of aggregation are well captured.

The fitted ICNC${}_{t_{\mathrm{0}}}$ exponent (β₀=0) is below 1, which is consistent with the above free power-law fitting, which includes only one controlling factor, ICNC${}_{t_{\mathrm{0}}}$. This strong, but nonlinear, positive contribution aligns with both the causal graph (Fig. 7) and SHAP analysis (Fig. 9), which also identified ICNC${}_{t_{\mathrm{0}}}$ as the dominant predictor. The temperature coefficient (β₁=0) indicates a positive sensitivity to temperature, consistent with laboratory evidence (Hosler et al., 1957) as well as the causal graph and SHAP findings, though the effect is modest. However, this coefficient reflects only the behavior within our observational temperature range (between −4.7 and −7.8 °C). Finally, the major size exponent (β₃=0) is relatively small, suggesting a weaker dependence on crystal size compared to ICNC${}_{t_{\mathrm{0}}}$ or temperature. This may reflect either the limited variation in crystal size across experiments, especially the colder cases (Fig. 2e) or the dominant influence of concentration effects under high-ICNC${}_{t_{\mathrm{0}}}$ seeded conditions.

4.4 Sensitivity to temperature and ICNC${}_{t_{\mathrm{0}}}$

We evaluated the sensitivity of predicted aggregation rates to temperature and ICNC${}_{t_{\mathrm{0}}}$ by comparing predictions from CatBoost (CatB) and the physical equation against observations (Fig. 10). Three prescribed ICNC${}_{t_{\mathrm{0}}}$ levels – 10¹, 10², and 10³ L⁻¹ – were selected to represent conditions characteristic of stratus clouds (Gultepe et al., 2001), deep convection (Heymsfield and Willis, 2014), and secondary ice production in convective cloud systems (Korolev et al., 2020). To enable a meaningful comparison between predictions and observations, we stratified the measurements into three corresponding ICNC_t1 intervals: (10⁰, 10¹), (10¹, 10²), and (10², 10³) L⁻¹. These intervals approximately align with the predicted ICNC${}_{t_{\mathrm{0}}}$ levels and allow an assessment of the models’ ability to capture the observed magnitudes and trends across representative regimes.

Both models reproduce the observed positive temperature dependence within the measurement domain (−7.8 to −4.7 °C; shaded region in Fig. 10) and exhibit physically plausible trends beyond. At all fixed ICNC${}_{t_{\mathrm{0}}}$, predicted aggregation rates increased systematically with temperature, broadly consistent with the observations. Observed rates generally clustered around the predicted values at their respective ICNC levels, supporting the proposed scaling.

CatBoost predictions were smooth and closely aligned with observations at high and intermediate ICNC levels, but at low ICNC they exhibited pronounced fluctuations at warmer temperatures (T> −6 °C) and a nearly flat response at colder temperatures (T< −6 °C; solid dark blue line in Fig. 10). This behavior reflects the high variability of observations in the low-ICNC regime. We quantify this variability using the coefficient of variation (CV), defined as $\mathrm{CV}=\mathit{\sigma }/\mathit{\mu }$, where σ and μ are the standard deviation and mean of R_agg within a given bin. Across the observed temperature range, CV values at low ICNC are typically 1.2–2.1, compared to 0.7–1.0 at intermediate ICNC and ≲ 0.6 at high ICNC. Such large CV values indicate that aggregation rates in this regime are dominated by strong fluctuations rather than a clear temperature signal. In combination with the imposed lower bound on the target variable during training, this lack of a resolvable signal causes the tree-based model to revert to the lower-bound baseline, resulting in an approximately constant prediction at low temperatures. The flat behavior therefore reflects the absence of a learnable temperature dependence in this regime, rather than a physically meaningful relationship. The fact that a localized dip near −5.3 °C appears consistently across all ICNC regimes, while being most pronounced at low ICNC and progressively weaker at higher ICNC, further supports this interpretation. This feature is more likely to reflect experiment-to-experiment variability and condition mixing near this temperature, together with the non-smooth, threshold-based nature of tree models, rather than a robust, generalizable physical transition.

In contrast, predictions from the physical equation were stable and monotonic across the full ICNC and temperature ranges, reproducing the general trends well, though slightly overestimating aggregation rates at the lowest ICNC and T< −7 °C. The increased variability in low-ICNC observations suggests that other factors, such as habit composition, size distribution, or turbulence, may exert a proportionally stronger influence when ICNC is low. From a modeling perspective, this implies that physically based parameterizations provide a more reliable baseline representation of aggregation under such conditions, with observational or machine-learning-based approaches best suited for diagnosing variability or informing uncertainty estimates rather than defining deterministic temperature dependence.

Although evaluation based on RMSE and R² over the full test set indicated slightly better performance of CatBoost compared to the physical equation, Fig. 10 suggests that the two approaches perform similarly in most regimes, with the physical equation outperforming CatBoost in some cases (e.g., low ICNC at T< −7 °C). This apparent discrepancy can be partly attributed to how the performance metrics are computed: RMSE and R² are dominated by data-rich regions, where CatBoost tends to follow the observed scatter more closely – especially when the observations themselves are subject to variability from binning into three ICNC_t1 intervals [(10⁰, 10¹), (10¹, 10²), and (10², 10³) L⁻¹]. In such regions, CatBoost can achieve lower residuals by adapting flexibly to local variations, whereas the physical equation prioritizes smooth, process-consistent behavior. In contrast, in sparsely sampled or more physically complex regimes, the physical equation yields more stable and monotonic trends, while CatBoost shows larger deviations. These complementary behaviors underscore the value of combining data-driven flexibility with physically constrained formulations.

Figure 10 Sensitivity of predicted and observed aggregation rates to temperature across ICNC levels. Predictions from the CatBoost model (CatB) and a physical equation are shown for three prescribed initial ice crystal number concentrations (ICNC${}_{t_{\mathrm{0}}}$) of 10¹ (solid), 10² (dashed), and 10³ L⁻¹ (dash-dotted); blue and red curves correspond to CatB and the physical equation, respectively. Observations are divided into three ICNC_t1 intervals: (10⁰, 10¹), (10¹, 10²), and (10², 10³) L⁻¹, and are plotted as crosses in cyan, orange, and red, respectively. The shaded region indicates the observed temperature range (−7.8 to −4.7 °C).

5 Conclusions

In this study, we present the in-situ quantification of ice aggregation rates in persistent supercooled stratiform clouds. This was made possible by a novel experimental design combining glaciogenic cloud seeding, which provided a controlled initial state and allowed us to determine the age of ice crystals from nucleation to observation, and IceDetectNet, a deep-learning algorithm capable of counting the number of monomers in individual ice crystals. This approach offered new insights into the microphysical and meteorological conditions that govern aggregation. The ICNCs observed during the experiments exceeded those typically found in naturally occurring stratiform clouds but were comparable to levels observed during secondary ice production events in convective clouds (Korolev et al., 2020). The combination of high ICNC and the absence of secondary ice production provided a rare opportunity to isolate the contribution of aggregation without the confounding influence of ice multiplication processes and to generate observational constraints relevant to natural clouds.

We found ICNC${}_{t_{\mathrm{0}}}$, temperature, major axis length, and aspect ratio as the primary controls on aggregation rate among all the factors we investigated, with ICNC${}_{t_{\mathrm{0}}}$ emerging as the dominant factor. Even though ICNC${}_{t_{\mathrm{0}}}$ has the strongest influence, a linear model based solely on ICNC${}_{t_{\mathrm{0}}}$ failed to reproduce the observed rates, showing the importance of other factors and the nonlinearity of their interactions. While ICNC${}_{t_{\mathrm{0}}}$, temperature, and major axis length contribute positively to aggregation, aspect ratio acts as a negative factor, and riming showed no detectable effect, showing that riming and aggregation are largely independent processes. The influence of temperature operates both directly and indirectly by modifying ice crystal shape, consistent with current microphysical schemes (e.g. Seifert and Beheng, 2006). EDR showed no significant correlation with the aggregation rate. This finding likely reflects the combined effects of narrow size distributions and the dominant roles of ICNC${}_{t_{\mathrm{0}}}$ and temperature in these seeded clouds. Within the EDR range sampled, aggregation appears to be at most weakly sensitive to EDR, and any residual dependence is small compared with the primary controlling factors.

We demonstrated that both machine learning models and the physically derived equation successfully reproduced the observed aggregation rates. CatBoost achieved the best performance in terms of RMSE and R² on the test dataset by capturing nonlinear interactions, whereas the physically based model proved more robust and stable in sensitivity tests, particularly in regimes with low ICNC${}_{t_{\mathrm{0}}}$, higher observational variability, and outside the observed temperature range. These differences highlight the complementary nature of data-driven and physically constrained approaches.

This study advances our understanding of the microphysical and environmental controls on ice aggregation. It emphasizes the central role of ICNC${}_{t_{\mathrm{0}}}$ and the influence of temperature-dependent ice crystal properties. The study also showed that the aggregation occurs within 5–10 min after the ice formed. However, the mechanism behind the observed sub-quadratic relationship between ICNC${}_{t_{\mathrm{0}}}$ and aggregation rate remains uncertain. One possible explanation that is consistent with our observations is that aggregation may involve smaller ice crystals. However, this remains hypothetical. Addressing these uncertainties requires experimental designs that can capture the intermediate stages of ice crystal growth and interaction. The single-point experimental setup limits our ability to capture the full evolution of microphysical processes and may obscure intermediate stages of aggregation. Future experiments would benefit from a Lagrangian observational setup to better resolve the complete aggregation pathway. This would allow us to refine the representation of aggregation in weather and climate models, ultimately improving predictions of precipitation and cloud radiative effects.

Appendix A: IceDetectNet-CLOUDLAB Uncertainties

To assess the uncertainty of IceDetectNet-CLOUDLAB, we compared its predicted monomer counts with manual annotations on a test set comprising 246 images of aggregates, totaling 359 labeled monomers. Discrepancies were calculated as the difference between predicted and annotated monomer counts. For example, if the hand label (HL) assigned two monomers while the prediction assigned one, the discrepancy was recorded as −1; if the prediction assigned three, the discrepancy was +1. Across all test images, the total discrepancy summed to 14, indicating a slight net overestimation of 14 monomers. To investigate the model's behavior across different aggregation levels, we grouped ice crystals by their hand-labelled monomer count (1 to 4) and computed the mean and standard deviation of the predicted counts within each group. Ideally, a perfect model would yield a mean equal to the true monomer count and a standard deviation of zero. The largest relative deviation occurred in the 2-monomer category, where the model overestimated the monomer count by 21.9 % (Table A1). Predictions for 3-monomer aggregates showed a slight underestimation. For the most frequent class 1-monomer ice, ice crystals – IceDetectNet-CLOUDLAB performed reliably, with a mean prediction of 1.10 and a low standard deviation of 0.35 (Table F1). These results indicate that the model is robust for monomer number counting but tends to slightly overestimate complexity in multi-monomer aggregates.

Table A1 Prediction uncertainty of IceDetectNet-CLOUDLAB. Comparison between predicted (Pred) and hand-labeled (HL) monomer counts on the test set. Discrepancy (Discr.) is defined as prediction minus ground truth. Discr. (%) is normalized by HL Count.

Appendix B: Correlation Between Aggregation Rate and EDR

We evaluated whether turbulence intensity, represented by EDR, influenced the aggregation rate across the 21 seeding experiments. No significant correlation was found in either the warmer (r=0, p=0) or colder ($r=-\mathrm{0}$, p=0) temperature regimes. Mean aggregation rates remained largely invariant across EDR bins, with substantial within-bin variability. HOLIMO imaged cloud particles within a three-dimensional volume of 11.76 cm³ at 20 Hz during seeding conditions (Fuchs et al., 2025; Ramelli et al., 2020, 2024), yielding an effective spatial resolution of approximately 1 m along the flight path.

Several factors may contribute to the absence of a detectable EDR signal: (1) the ice crystal size distribution in seeding experiments was relatively narrow and often dominated by specific ice habits, unlike the broader and more complex size distributions typical of natural mixed-phase clouds. This simpler distribution composition may reduce the sensitivity of aggregation processes to turbulence; and (2) the dominant influence of ICNC and temperature on aggregation rates could obscure turbulence effects, particularly in short-lived seeded clouds. A more straightforward interpretation is that aggregation is only weakly sensitive to EDR within the range of dissipation rates sampled here, and any remaining dependence is small relative to the dominant effects of ICNC${}_{t_{\mathrm{0}}}$ and temperature.

Figure B1 Correlations between aggregation rate and EDR grouped by temperature. Same as Fig. 2c, but with EDR as the x axis.

Figure B2 Distribution of EDR grouped by temperature group. Top: Probability density functions of EDR for all seeding experiments, grouped by temperature: colder (T≤ −7 °C, blue) and warmer (T> −7 °C, red). Bottom: Mean EDR values for individual experiments, shown as red circles (warmer) and blue squares (colder).

Appendix C: Machine Learning Aggregation prediction Models

We used a total of 10 supervised regression algorithms. A general introduction to these machine learning methods can be found in Géron (2022). Below we summarize the basic principles of each model:

• Linear Regression fits a linear function by minimizing the squared difference between predicted and observed values. It assumes independent, additive relationships among input features and cannot capture nonlinear interactions.

• Ridge Regression extends linear regression by adding L₂ regularization, which penalizes large coefficients and reduces overfitting, particularly in the presence of small sample sizes.

• Bayesian Ridge Regression introduces a probabilistic prior over the regression coefficients and estimates them using Bayesian inference. This yields both regularized predictions and uncertainty estimates, while still assuming linearity.

• K-Nearest Neighbors is a nonparametric method that predicts output values by averaging the labels of the k training samples closest to the test point. Similarity is typically measured by Euclidean distance in the standardized feature space. Although simple and interpretable, the method becomes less effective in high-dimensional spaces.

• Decision Tree Regression recursively partitions the input space by choosing feature-value thresholds that minimize the prediction error (typically mean squared error) at each split. The resulting model is a set of hierarchical “if–then” rules – for example, “if ICNC${}_{t_{\mathrm{0}}}$ > 800 L⁻¹ and temperature < −6 °C, then…”. At each node, the best split is selected without looking ahead. While decision trees can capture nonlinear relationships and feature interactions, they tend to overfit if not regularized.

• Adaptive Boosting constructs an ensemble of weak learners – models that perform only slightly better than random guessing – by sequentially reweighting the training samples. After each iteration, higher weights are assigned to mispredicted samples, forcing the next learner to focus more on difficult cases. The final output is a weighted sum of all learners. While this can substantially reduce bias, AdaBoost can be sensitive to noise and outliers.

• Gradient Boosting also builds an ensemble of decision trees, but instead of adjusting sample weights, it fits each new tree to the residuals (i.e., errors) of the combined ensemble so far. This stage-wise additive approach allows the model to incrementally improve predictions and capture complex nonlinear dependencies.

• Light Gradient Boosting Machine (Ke et al., 2017) is an optimized gradient boosting implementation. It uses histogram-based feature binning, where continuous input features are discretized into fixed-width bins to accelerate training and reduce memory usage. Unlike traditional level-wise tree growth, LightGBM grows trees leaf-wise by expanding the leaf with the highest loss, leading to deeper, more specialized trees. These optimizations make LightGBM highly efficient on large datasets.

• Extreme Gradient Boosting enhances standard gradient boosting with L₁ and L₂ regularization to penalize model complexity and prevent overfitting. It also supports parallelized training and handles missing values natively during tree construction. These improvements make it robust and scalable for structured data tasks.

• Extremely Randomized Trees is an ensemble of decision trees where both the features and split thresholds are selected randomly at each node, rather than chosen based on an optimal impurity measure. This high degree of randomness reduces variance and helps avoid overfitting, especially when the data contain noise.

Appendix D: Correlation Between Aggregation Rate and ICNC

Figure D1 ICNC${}_{t_{\mathrm{0}}}$ – aggregation rate relationships across all experiments. Same as Fig. 3, but shown individually for all 21 experiments. Red and blue markers represent warmer (T> −7 °C) and colder (T≤ −7 °C) cases, respectively.

Appendix E: Examples of ice crystals

Figure E1 A randomly selected sample of ice crystal images observed by HOLIMO during (a) SM054 and (b) SM069.

Appendix F: Experiments overview

Table F1 Summary of all seeding experiments. For each experiment (Exp.) (where “SM” denotes seeding mission, consistent with Fuchs et al., 2025; Miller et al., 2025), the table lists the number of 1 s data points (#Data), advection time between seeding and sampling (Time, in s), wind speed at seeding level (Wind, in m s⁻¹), temperature (Temp, in (°C)), mean ice crystal number concentration (ICNC, in cm⁻³), background liquid water content (LWC, in mg m⁻³), cloud droplet number concentration (CDNC, in cm⁻³), and mean eddy dissipation rate (EDR, in 10⁻³ m² s⁻³). ICNC, LWC, CDNC, and EDR reported values are means ± standard deviations.