Latent heat feedbacks and the self-lofting of seeded ice plumes: insights from bin microphysics simulations

The interaction between environmental dynamics and microphysical processes governs the efficacy of glaciogenic cloud seeding, yet the extent to which initial vertical wind conditions dictate the evolution of seeded ice plumes remains poorly constrained. We investigate the dynamical and microphysical life cycle of ice plumes in supercooled stratiform clouds using idealized Large-Eddy Simulations (LES) coupled with the bin microphysics scheme SCALE-AMPS. Simulations of a targeted CLOUDLAB seeding experiment are constrained and validated by in situ observations. An ensemble of 20 simulations initialized with varying vertical wind velocities reveals a fundamental transition: while the initial trajectory is kinematically governed by the environmental vertical wind at cloud seeding, the long-term evolution is dominated by internal thermodynamic feedbacks. We identify a “self-lofting” mechanism wherein buoyancy generated by latent heat release from rapid ice growth overcomes initial subsidence, causing even downdraft-seeded plumes to eventually ascend. This thermodynamic response creates a structural trade-off: plumes initiated in updrafts are terminated by the cloud-top inversion, whereas those in downdrafts experience delayed ascent, resulting in a lager vertical dispersion. This expanded vertical extent compensates for the lower mean altitude of downdraft plumes, ensuring their geometrical detectability by downstream sampling in in observations. Microphysically, downdraft plumes undergo transient sublimation near cloud base but recover following buoyancy-driven re-ascent. These findings demonstrate that glaciogenic seeding is dynamically robust in low stratus clouds, as the seeded ice plume acts as an active thermodynamic agent capable of sustaining its residence time in the mixed-phase layer independent of the initial vertical wind state.

1 Introduction

Ice microphysical processes govern the phase partitioning, persistence, and radiative properties of mixed-phase clouds, thereby exerting a critical influence on Earth’s hydrological cycle (Pruppacher et al., 1998; Korolev et al., 2017). These processes are intrinsically coupled to atmospheric dynamics: the interplay between vertical air motion, turbulence, and local supersaturation dictates ice crystal growth modes and cloud longevity (Korolev et al., 2003; Shupe et al., 2008). Vertical air motions are widely recognized as a primary driver of this evolution (Bühl et al., 2019; Shupe et al., 2008). Updrafts promote condensation and maintain liquid water by counteracting the Wegener–Bergeron–Findeisen (WBF) process – a mechanism wherein ice crystals grow at the expense of evaporating supercooled liquid droplets. In contrast, downdrafts and subsidence favor glaciation and liquid water depletion (Rauber and Tokay, 1991; Korolev and Field, 2008; Bühl et al., 2019).

While vertical velocity fluctuations modulate ice nucleation efficiency and particle trajectories (Pinsky and Khain, 2002; Grabowski and Abade, 2017; Abade and Albuquerque, 2024), subsequent ice growth is accompanied by latent heat release. This thermodynamic response alters buoyancy and feeds back onto vertical motions, thereby further coupling microphysical and dynamical evolution (Morrison et al., 2012). Following the initial ice initiation – the so-called “latent heat kick” from droplet freezing – the ice phase undergoes rapid diffusional growth. During this stage, the latent heating produced by depositional growth significantly outweighs the additional weight of the ice phase (condensate loading), effectively invigorating the cloud plume. While some studies have explored the resilience of cloud microphysics to environmental changes (Grabowski and Morrison, 2021), the specific roles of the initial freezing “kick” versus the subsequent feedbacks driven by the WBF mechanism and adiabatic cooling in governing plume dynamics remain under-examined. A fundamental question, therefore, arises: to what extent is the evolution of an ice plume determined by the initial vertical wind conditions at the time of nucleation, as opposed to being sustained and amplified by the seeding-induced microphysical – thermodynamic feedbacks – including both the initial latent heat “kick” from droplet freezing and the subsequent latent heat release associated with the transition from liquid and vapor to the ice phase?

Shallow mixed-phase clouds, such as stratocumulus and altocumulus, serve as an ideal natural laboratory for investigating these interactions. Compared to deep convection, these clouds are dynamically simpler, yet they exist in a thermodynamically metastable state where microphysical evolution is highly sensitive to the partitioning between liquid and ice phases (Ansmann et al., 2009; Westbrook and Illingworth, 2013). Glaciogenic cloud seeding offers a powerful framework for interrogating these systems by introducing a controlled source of ice crystals into supercooled clouds at a known time and location. Historically, glaciogenic seeding has been extensively studied to assess its potential for precipitation enhancement, with major field campaigns such as SNOWIE providing comprehensive evidence of cloud-seeding-induced snowfall over mountainous terrain (Tessendorf et al., 2019; Friedrich et al., 2021). While these meso-scale experiments have successfully demonstrated the macroscopic impact of seeding, the fundamental microphysical-dynamical interactions at the plume scale remain challenging to isolate. The CLOUDLAB (Henneberger et al., 2023) campaign employing uncrewed aerial vehicles (UAVs) and in situ instrumentation has enabled detailed characterization of seeded ice plumes, revealing rapid ice growth, efficient liquid water depletion, and enhanced aggregation (Miller et al., 2024, 2025; Ramelli et al., 2024; Fuchs et al., 2025; Zhang et al., 2026).

Nevertheless, such measurements are inherently Eulerian, providing snapshots of plume properties at fixed locations rather than a continuous record of their evolution. As a result, distinguishing whether the vertical transport of an ice plume is driven by the ambient vertical velocity into which it was born or by the buoyancy generated by its own growth remains difficult to establish from observations alone. Large-eddy simulations (LES) coupled with detailed bin microphysics schemes provide the necessary complementary tool to bridge this gap. Models such as the Scalable Computing for Advanced Library and Environment–Advanced Microphysics Prediction System (SCALE-AMPS; Ong et al., 2022) explicitly represent the size-, shape-, and habit-dependent evolution of ice particles while capturing the feedback between phase changes and air motion (Hashino and Tripoli, 2007).

In this study, we use the SCALE–AMPS model to simulate a targeted seeding event from the CLOUDLAB campaign, utilizing in situ observations to constrain the initial conditions and validate the microphysical evolution. The paper is organized as follows. We first describe the experimental setup and the in situ dataset (Sect. 2), followed by the numerical configuration of the SCALE–AMPS model (Sect. 3). The model performance is then evaluated against observations by comparing key parameters such as Ice Crystal Number Concentration (ICNC), Liquid Water Content (LWC), and particle size distributions (Sect. 4.1). Subsequently, we present the core sensitivity analysis, isolating the impact of initial vertical wind conditions on plume evolution and quantifying the role of latent heat feedbacks (Sect. 4.2 and 4.4). Finally, the main findings and their implications for future mixed-phase cloud studies are summarized (Sect. 5).

2 Observational Setup and Data

This study focuses on a single glaciogenic seeding experiment from the CLOUDLAB campaign in persistent wintertime supercooled stratus over the Swiss Plateau near Eriswil, Switzerland (47°04^′14^′′ N, 7°52^′22^′′ E; 920 m a.s.l.) (Scherrer and Appenzeller, 2014; Henneberger et al., 2023; Miller et al., 2024, 2025; Ramelli et al., 2024; Fuchs et al., 2025). The experiment was conducted on 24 January 2023 at 18:08:57 UTC in liquid-dominated, quasi-stationary cloud conditions with temperatures close to −5 °C at seeding height. Glaciogenic seeding was performed at an altitude of 1348 m, 2068 m upwind of the measurement site, under mean wind directions between 75 and 77°. As the seeded ice plume advected over a fixed ground-based site, its microphysical changes were continuously sampled by an in situ platform, providing time-resolved observations used both to constrain the model setup and to evaluate the simulated ice growth and interaction processes. The sampling altitude during plume observation was generally consistent with the seeding height, with vertical fluctuations of up to about ± 30 m.

In situ measurements were obtained using a tethered balloon system carrying the HOLographic Imager for Microscopic Observations (HOLIMO; Ramelli et al., 2020), a digital in-line holography instrument with effective spatial resolutions of 6 µm for cloud droplets and 25 µm for ice crystals. Holograms were reconstructed, and detected particles were initially classified using a convolutional neural network (Touloupas et al., 2020). Ice crystals classified by this procedure were subsequently analyzed using IceDetectNet-CLOUDLAB (Zhang et al., 2026), a CLOUDLAB-specific fine-tuning of IceDetectNet (Zhang et al., 2024). The uncertainty in cloud droplet number concentration is approximately 5 %, while uncertainties in ice crystal number concentration range from 5 %–10 % for larger crystals and increase to about 15 % for smaller particles.

The numerical simulation is configured to reproduce this experiment using observationally constrained initial and environmental conditions as follows: The model is initialized with a liquid water path of 107.5 g m⁻² measured from HATPRO Microwave Radiometer. Thermodynamic profiles for the model setup are derived from radiosonde launches at 15:44, 20:31, and 20:48 UTC on 24 January 2023, which span the seeding and plume observation period. Model performance is evaluated against the observed ice size distribution, aspect ratio, ICNC, CDNC, LWC, and the ice–ice collision rate inferred from observed ice microphysical properties and deep learning techniques (see Zhang et al., 2026 for detailed methodology), using consistent definitions and averaging procedures for model output and observations.

3 Method

3.1 Model

The SCALE-AMPS model (Nishizawa et al., 2015; Sato et al., 2015; Ong et al., 2022) is used in this study. SCALE is a large-eddy simulation (LES) atmospheric model that solves the primitive equations using the finite-volume formulation on a structured Arakawa C-staggered grid. Prognostic variables are advanced in time using the fourth-order Runge–Kutta integration scheme. Subgrid-scale turbulence is represented using a Smagorinsky–Lilly–type scheme (Smagorinsky, 1963; Lilly, 1962), and radiative transfer is computed with the MSTRNX radiation code (Sekiguchi and Nakajima, 2008). The current subgrid-scale scheme does not account for local supersaturation fluctuations (e.g., Grabowski and Abade, 2017). This idealization implies that microphysical processes, such as ice nucleation and droplet activation, are computed based on grid-mean thermodynamic states, a common simplification in current LES frameworks.

AMPS is a habit-preserving bin microphysics scheme. Detailed descriptions of AMPS are provided in Hashino and Tripoli (2007, 2008, 2011), and its implementation within the SCALE framework is described by Ong et al. (2022). Here, only a brief overview is given. Both liquid and ice particle size spectra are discretized into 40 bins. For each liquid bin, tracer variables include the liquid water mixing ratio (Q_liq), Cloud Droplet Number Concentration (CDNC), aerosol mixing ratio (Q_ae), and the mean ratio of soluble mass to total aerosol mass (Q_solv). For each ice bin, the corresponding variables are the ice mixing ratio (Q_ice), ICNC, Q_ae, and Q_solv.

To keep track of the evolution of the shape of ice particles, three additional tracer variables corresponding to the a-axis and c-axis lengths are stored in each ice bin. These variables allow the model to keep track of the time evolution of habit changes, as ice particles grow into different shapes (e.g., planar or columnar) depending on ambient temperature. Aggregation (Q_agg) and riming (Q_rim) mixing ratios are also predicted, because collisions with ice particles and liquid droplets alter particle mass compositions and shapes. AMPS assumes that aggregation and riming processes cause ice particles to asymptotically approach an elliptical shape. A third mass-component prognostic variable, the crystal mass mixing ratio (Q_cry), is included to track the history of vapor deposition growth. Vapor deposition growth increases only Q_cry, while Q_agg and Q_rim remain unchanged. For pristine ice particles that grow exclusively by vapor deposition, Q_cry=Q_ice within a bin.

A cubic polynomial is used to represent the subgrid-scale distribution of number density in both the liquid and ice spectra. During diffusional or collisional growth, particle mass changes cause the lower and upper boundaries of a bin to shift. The subgrid distribution in each shifted bin is recalculated, and the prognostic variables in the original bins are updated by integrating the corresponding quantities over the overlapping regions between the shifted and original bins. The ratios Q_cry $/$ Q_ice, Q_agg $/$ Q_ice and Q_rim $/$ Q_ice are assumed to remain constant within the subgrid distribution. This is a strategy called implicit mass sorting assumption (Hashino and Tripoli, 2007). When two ice particles collide, the post-collision axis lengths are computed using the volumetric means.

The relative growth rates of the a- and c-axis lengths are determined by an inherent growth-rate function that depends solely on ambient temperature (Chen and Lamb, 1994). AMPS uses a single spectrum for all types of ice particles. Ice particle types can be diagnosed within each bin based on the crystal, aggregation and riming mass fractions (Q_cry $/$ Q_ice, Q_agg $/$ Q_ice, and Q_rim $/$ Q_ice). In AMPS, ice particles are classified as one of the following types: pristine crystals, rimed crystals, graupel, aggregates, or rimed aggregates (Hashino and Tripoli, 2007). For example, ice particles are diagnosed as aggregates when the aggregation mass fraction exceeds 10 % and the riming mass fraction is less than 10 %.

Cloud condensation nuclei (CCN) are treated as a spatially uniform background variable throughout this study and are assumed to follow a lognormal size distribution. The distribution is characterized by a total number concentration of 317 cm⁻³ (corresponding to the observed mean cloud droplet number concentration), a geometric mean radius of 0.052 µm, and a geometric standard deviation of 0.71. The CCN are assumed to consist entirely of ammonium bisulfate (NH₄HSO₄). The distribution is discretized into 10 bins in the computation of CCN activation. CCN in bins with critical supersaturations lower than the ambient supersaturation are activated and transferred to the liquid water bins, with the resulting droplet size and bin assignment determined by the microphysics time step and the initial droplet growth rate. Activated CCN are immediately replenished to maintain the prescribed background distribution.

Figure 1 Vertical profiles of (a) temperature, (b) wind speed, and (c) liquid water content (LWC) used as the initial conditions for the simulations. The mean wind is directed in the positive y direction. Panel (d) shows the imposed large-scale subsidence W_sub. The gray shaded region denotes the cloud layer.

Figure 2 1 min running average of the vertical wind speed in the grid cells where ice-nucleating particles (INPs) are released. The dots indicate the specific vertical wind speeds at the 20 distinct initiation times for the ensemble members. (b) schematic illustration of the computational domain after the spin-up period. Cloud seeding is performed between the coordinates (950, 10, 310) and (1050, 10, 310) m.

3.2 Numerical setup

The computational domain spans 2 and 4 km in the x and y directions, respectively, with the model top at 1 km. The grid spacing is uniformly 20 m in the x, y, and z directions. Periodic boundary conditions are applied in the horizontal directions, and Rayleigh damping is imposed above 800 m. Complex orography over the Swiss Plateau at the field campaign site is not represented; consequently, the domain is flat. This idealization may lead to an underestimation of turbulent mixing associated with real terrain. The cloud layer extends from approximately 220 to 600 m in altitude. Accordingly, the relative humidity is prescribed as 100 % within this layer and varies linearly below and above the cloud, with boundary conditions derived from radiosonde observations. The initial vertical profiles of temperature and liquid water content (LWC; Fig. 1) are derived from the linearized relative humidity profile, observed pressure, and the assumption of a vertically uniform liquid potential temperature below cloud top. The initial wind profile is obtained by smoothing the radiosonde observations, with the mean flow oriented in the y direction. In all simulations, the wind speed is nudged toward the initial profile with a relaxation time scale of 1 h. Model output is saved at 12 s intervals.

A constant large-scale subsidence of 0.00952 m s⁻¹ (Fig. 1d) is imposed above the inversion layer at 600 m above ground to maintain a steady boundary layer depth and inversion strength throughout the simulations, as the simulations are idealized and do not include realistic boundary conditions. The subsidence is assumed to be zero at the surface and is linearly interpolated at altitudes between the inversion layer and the surface. Its effect on temperature and water vapor mixing ratio is represented by adding a source term $-w\partial \mathit{\varphi }/\partial z$ in the energy and moisture equations, where w is the subsidence and ϕ denotes temperature or water vapor mixing ratio, to their respective prognostic tendencies.

Observational constraints on the physical and chemical properties of background CCN are unavailable. We therefore assume that the CCN consist entirely of ammonium bisulfate, with a mean particle radius of 0.052 µm and standard deviation of the CCN size distribution of 0.071. The observed mean cloud droplet number concentration is 317 cm⁻³; accordingly, the background CCN number concentration is prescribed as 317 cm⁻³.

The first hour of integration is treated as a spin-up period to allow the model to adjust dynamically and thermodynamically to the imposed initial and boundary conditions. A spin-up time of 1 h is sufficient for turbulent flows in the boundary layer to fully develop and for the mean liquid water path and CDNC to approach quasi-steady values (not shown). After the spin-up period, a series of updrafts and downdrafts emerges within the boundary layer (see the schematic illustration of the idealized simulation setup in Fig. 2). These motions are generated by thermal convection in the presence of wind shear and surface friction. They propagate downstream and evolve with time, with existing updrafts and downdrafts dissipating and new ones forming. Their vertical extent can reach the cloud top.

Ice-nucleating particles (INPs) are released within grid cells spanning the coordinates (950, 10, 310) to (1050, 10, 310) m, where the three indices denote the x, y, and z coordinates, respectively. Because the INP release rate is not constrained by observations, we prescribe an injection duration of 1 min and a spatially uniform release rate of 10 cm⁻³ s⁻¹, such that the simulated ice crystal number concentration (ICNC) matches the observed ICNC in the reference simulation described later. The injection altitude is chosen to match the actual seeding altitude, at which the ambient air temperature upon seeding is −5 °C. Ice particles form from the injected INPs and are subsequently advected downstream by the mean wind.

The freezing efficiency of INPs is parameterized following Omanovic et al. (2024), with the frozen fraction defined as

$$\begin{array}{}\text{(1)}& F_{\mathrm{ice}}={\displaystyle \frac{-b}{\mathrm{1}+\exp (-k(T-T_{\mathrm{0}}\left)\right)}}+b\end{array}$$

where F_ice is the frozen fraction, T is temperature (K), b=0.97, k=0.88, and T₀=263.95 K. New ice particles form when F_iceN_INP, where N_INP is the total INP number concentration in a grid cell, exceeds the existing ICNC. All INPs are assumed to have a uniform diameter of 1.0 µm and are treated as passive tracers. Consequently, they do not interact dynamically with the flow, and their chemical composition and mass are not explicitly represented.

After the spin-up period, the 1 min running mean vertical wind speed within the seeding grid cells fluctuates between −0.6 and 0.6 m s⁻¹ (left panel of Fig. 2). To investigate the sensitivity of the ice plume evolution to vertical wind speed, we perform 20 seeding simulations, with the initiation time of INP release staggered at 1 min intervals. The simulation in which INPs are released 9 min after the spin-up period is selected as the reference case, and its results are compared with observational data to evaluate the validity of the model and the numerical setup. The 20 simulations are categorized into three groups: (1) updraft, (2) no-wind, and (3) downdraft, based on the mean vertical wind speed at the seeding location greater than 0.2 m s⁻¹, between −0.2 and 0.2 m s⁻¹, and less than −0.2 m s⁻¹, respectively. This classification yields 7 updraft, 6 no-wind, and 7 downdraft simulations. For the sensitivity analysis, the ice plume is defined as the set of grid cells in which the ICNC exceeds the 70th percentile, corresponding to ICNC values of approximately >100 L⁻¹). This definition is adopted to more effectively track the core of the ice plume as it disperses over a larger region during the later stages of plume evolution, although it differs from the observational definition, in which the plume is identified using an ICNC threshold exceeding 10 L⁻¹. The conclusions drawn from the analysis are not sensitive to this difference, except that the simulated vertical plume spread and standard deviations become larger.

To further isolate the physical mechanisms driving the plume evolution and to address the relative importance of the initial latent heat “kick” versus the subsequent WBF-driven feedbacks, two additional sets of 20-member ensemble simulations were conducted. These 40 simulations use the same initiation times as the original seeding ensemble to ensure a direct comparison under identical dynamical backgrounds: (1) WBF-off seeding case (L_s=L_c): In these simulations, the latent heat of sublimation (L_s) is set equal to the latent heat of condensation (L_c) in the model's thermodynamic equations. This modification effectively eliminates the net latent heat gain ($L_{\mathrm{f}}=L_{\mathrm{s}}-L_{\mathrm{c}}$) otherwise produced when ice crystals grow at the expense of cloud droplets (the WBF process). Furthermore, it ensures that any depositional growth driven by adiabatic cooling during the plume's ascent only releases heat equivalent to L_c, thereby removing the extra energetic contribution of freezing (L_f). This setup allows the initial latent heat release from the freezing of cloud droplets (the “initial kick”) to occur normally while suppressing the enhanced thermodynamic feedback during subsequent ice growth; (2) No-ice seeding case: In this baseline ensemble, INPs are released into the domain as passive tracers, but the ice nucleation process is disabled. This case serves as a reference to represent the plume's evolution governed solely by the background turbulent flow, in the absence of any microphysical–thermodynamic feedbacks. By comparing these two sensitivity ensembles with the original seeding simulations, we can quantitatively decompose the total buoyancy gain into contributions from the initial freezing event and the subsequent growth-driven intensification supported by both WBF and adiabatic cooling.

The mean wind speed in the y direction exhibits weaker variability than the initial profile by the end of the spin-up period, as the boundary layer becomes well mixed by convection and turbulence. The wind speed decreases only slightly from approximately 5.1 m s⁻¹ near the surface to 4.8 m s⁻¹ near the cloud top. Consequently, the ice plume in all simulations advects downstream at an approximately constant speed in the y direction (Fig. 3). In the reference simulation, the ice plume travels 2 km in 444 s. The time required for the ice plume to travel 4 km varies between 804 and 852 s among the simulations.

Figure 3 Time series of the downstream distance (y position) advected by the ice plume for all simulations.

Figure 4 Simulated and observed (a) ICNC, (b) CDNC, (c) LWC, and (d) collision rate at 2 km downstream of the seeding location. The simulation results are sampled at the height between 400 and 450 m above ground, while the observational data are taken at 400 m above ground.

4 Results

4.1 Model validation

In this section, the validity of the numerical setup is assessed by comparing results from the reference simulation with in situ observations obtained on 24 January 2023. Specifically, comparisons are made for ICNC, CDNC, LWC, ice crystal aspect ratio, and size distribution. The observational data used for comparison were collected using the HOLIMO instrument at an altitude of approximately 400 m above ground level (a.g.l.) and approximately 2 km downstream of the seeding drone (the balloon carrying the instrument was advected by the wind). The data were collected starting at 18:16:25 UTC, when the instrument first detected ice crystals. Each data point represents a time snapshot with an ice crystal present. The means, medians, and standard deviations are computed from this time series.

In the simulation, these quantities are sampled from grid cells located 2 km downstream of the seeding region and between 400 and 450 m a.g.l., when the ice crystal plume passes through these grid cells and the ICNC exceeds 10 L⁻¹. This sampling strategy is designed to match the location and conditions of the observational measurements. The means, medians, and standard deviations are computed in the same manner as for the observational data.

Figure 5 Simulated and observed (a) ice count and (b) aspect-ratio distributions at 2 km downstream of the seeding location.

Figure 6 IWC and LWC (contour lines) at x=2000 m at 4, 8, and 12 min after seeding begins.

The bin resolution in AMPS differs from that of the instrument. In AMPS, the size distribution is defined as a function of ice crystal mass and follows a logarithmic mass grid (Hashino and Tripoli, 2007). To enable direct comparison of the mean ice crystal aspect ratio and number concentration in each bin, the bin boundaries are converted from mass to equivalent spherical diameter using the volume relation of a sphere, assuming an ice crystal density of 0.916 g cm⁻³. The observed aspect ratio and size distribution data are then re-binned accordingly.

The simulated mean values of ICNC, CDNC, LWC, and collision rate agree reasonably well with the observed means (Fig. 4), considering that the setup is idealized. Although the simulated median CDNC and LWC are approximately 70 cm⁻³ and 0.05 g m⁻³ lower than the observed medians, owing to a higher simulated median ICNC, the interquartile range (25th–75th percentiles) of the simulated LWC lies within the range of the observed LWC. The observational data exhibit greater variability than the simulation. This discrepancy may be attributed to three factors: (1) small-scale inhomogeneity that is not fully resolved by the grid spacing used in the simulation; (2) the idealized setup, which employs smoothed initial conditions and neglects the complex terrain; and (3) the use of a deterministic ice nucleation scheme, which may underestimate the natural stochastic variability of ice initiation.

In the simulation, ice crystals grow predominantly by vapor deposition. Consequently, the simulated size distribution is narrower than that observed (Fig. 5). Most simulated ice crystals have sizes between 173 and 320 µm, whereas the observed size distribution peaks between 236 and 320 µm. This difference is also likely due to the absence of turbulent effects that would otherwise broaden the ice crystal size distribution (Grabowski and Abade, 2017), as well as insufficient spatial inhomogeneity.

For ice crystals with sizes between 173 and 320 µm, the simulated mean aspect ratio ranges from 4 to 5, in good agreement with observed values (Fig. 5). However, the simulated mean aspect ratio of ice crystals smaller than 173 µm reaches approximately 7. This bias arises because some ice crystals form at colder temperatures after the ice plume ascends to higher altitudes. Over this temperature range in the boundary layer, the aspect ratio (c-axis $/$ a-axis) increases more rapidly at colder temperatures (Chen and Lamb, 1994), leading to larger aspect ratios for ice crystals formed later in the plume evolution.

The ice plume initially appears as a localized feature, with a horizontal extent of approximately 300 m and a top height near 400 m 4 min after INPs are released into the computational domain (Fig. 6; see also Appendix A for 3D visualizations of the ice plume from the model simulations). After traveling 2 km downstream, the plume top reaches a height of about 500 m above ground, while its horizontal extent remains nearly unchanged. The vertical extent of the ice plume continues to increase with time, and it eventually spans the entire cloud depth at 3 km downstream.

Overall, although the simulated ice plume exhibits greater spatial homogeneity and potentially weaker turbulent effects than observed, the model is able to reproduce the key mean characteristics of the ice plume. Good agreement is achieved between the reference simulation and observations in terms of mean ICNC, CDNC, LWC, collision rate, ice crystal size, and aspect ratio. This numerical configuration therefore provides a foundation for the subsequent sensitivity experiments investigating the influence of vertical velocity at the seeding location.

4.2 Latent heat feedbacks and vertical plume dispersion

The vertical trajectory of the seeded ice plume is initially kinematically governed by the environmental vertical wind field at the point of injection. As illustrated in Fig. 7a, plumes injected into updraft regions undergo rapid ascent, reaching the cloud-top height (z≈600 m above ground) within 6 to 10 min. Conversely, plumes seeded into subsiding air parcels (w_init<0) experience an initial loss of altitude, descending by up to 150 m within the first 300 s. However, a fundamental dynamic transition is observed as the plume evolves: the trajectory decouples from the initial environmental forcing, and all plumes eventually establish a coherent net ascent. Even simulations initialized in strong downdrafts reverse their vertical motion, initiating a secondary ascent phase after approximately 300 to 500 s.

We hypothesize that this unconditional lofting is driven by internal microphysical-thermodynamic feedbacks, specifically the buoyancy generated by the latent heat of phase changes. The rapid growth of seeded ice crystals-primarily through vapor deposition in the supersaturated mixed-phase environment-releases a significant amount of latent heat. This exothermic process induces a localized positive temperature perturbation within the plume relative to the environment, thereby generating buoyancy that opposes and eventually overcomes the initial downdraft momentum.

Figure 7 Temporal evolution of ice plume dynamics and thermodynamic feedbacks. (a) Mean altitude of the ice plume center of mass. The observational sampling time corresponds to t≈464 s. (b) Evolution of the mean vertical wind speed (W) within the ice plume. (c) The relationship between the seeding-induced vertical velocity perturbation (ΔW) and the accumulative latent heat release. The strong positive correlation confirms that latent heating drives the secondary ascent.

Figure 8 Boxplots of the vertical extent (altitude spread) of the ice plumes sampled at downstream distances of (a) y=2, (b) y=3, and (c) y=4 km. In each boxplot, the central horizontal line and the white circle represent the median and the mean values, respectively. The top and bottom edges of the box indicate the 75th and 25th percentiles, while the whiskers extend to the 5th and 95th percentiles. Simulations are categorized by the vertical wind conditions at the injection time: Updraft (w_init>0.2 m s⁻¹), No vertical wind ($-\mathrm{0.2}\le w_{\mathrm{init}}\le \mathrm{0.2}$ m s⁻¹), and Downdraft ($w_{\mathrm{init}}<-\mathrm{0.2}$ m s⁻¹).

This mechanism is quantitatively substantiated by two lines of evidence. First, the mean vertical velocity (W) within the plume converges towards positive values across the entire ensemble, independent of the initial kinematic state (Fig. 7b). Second, and more critically, there is a robust functional relationship between the seeding-induced vertical velocity perturbation (ΔW) and the cumulative latent heat release (Q−Q_{no ice}), as shown in Fig. 7c. The cumulative latent heat is calculated as the time-integral of the net phase-change energy within the plume cells (defined as grid cells where ICNC exceeds the 70th percentile, corresponding to ICNC values of approximately 100 L⁻¹ or greater). This term specifically represents the excess heating relative to the non-seeded baseline (Q−Q_{no ice}), thereby isolating the latent heat of freezing and the net heating from the WBF process while excluding the background latent heat of condensation associated with the ambient cloud. The increase of ΔW with accumulated heat confirms that the secondary updraft is not a stochastic feature of the turbulence but a deterministic response to the latent heat released upon phase changes.

The interaction between this buoyancy-driven ascent and the thermodynamic stability at the cloud top (the capping inversion) fundamentally alters the vertical morphology of the plume. Plumes in updrafts possess high vertical momentum and reach the inversion layer rapidly (short residence time in the boundary layer). Upon impact, their ascent is stopped, causing the plume to spread horizontally (see Appendix A). In contrast, plumes in downdrafts exhibit a delayed ascent and longer residence time within the mixed-phase layer. This extended transit allows for greater vertical dispersion before reaching the ceiling. Consequently, as evidenced in Fig. 8, downdraft-seeded cases cause a significantly larger vertical spread (up to ∼250 m) at downstream distances compared to the vertically constrained updraft cases.

These dynamical differences have critical implications for the observational verification of the seeding signal. The in-situ measurement sampling was conducted at a downstream distance corresponding to an advection time of t≈464 s, at an altitude proximate to the seeding level. At this temporal snapshot (Fig. 8a), a significant vertical bifurcation is evident: while the centroids of the updraft-seeded plumes have ascended well above the sampling height (exceeding 500 m above ground), the downdraft-seeded plumes remain low, with centroids as low as ∼150 m. However, the enhanced vertical dispersion in the downdraft regime plays a critical compensatory role. These plumes occupy a much deeper vertical column, with a vertical extent reaching up to ∼250 m (Fig. 8a). This expanded vertical geometry ensures that even when the plume center of mass is below the flight level, the upper parts of the plume likely intersect the sampling height. Furthermore, this structural characteristic persists further downstream, as evidenced by the consistent spread statistics at y=3 and y=4 km (Fig. 8b, c). Consequently, despite the substantial variability in initial vertical trajectories, the seeding signature remains robustly detectable, as the thermodynamic feedbacks effectively broaden the spatial probability of interception.

Figure 9 Time evolution of the mean altitude of the ice plume in all simulations, with INP injection initiated at 1 min intervals from 0 to 1140 s (19 min). Solid curves represent the Normal configuration; dashed curves show results from simulations in which latent heating associated with the WBF process is disabled by setting the latent heat of sublimation equal to that of evaporation; dash-dotted curves indicate the mean altitude of INPs in simulations without ice nucleation.

4.3 Decomposing thermodynamic drivers: delayed feedback and the role of WBF

To disentangle the distinct thermodynamic drivers of the plume evolution, we compare the ice plume trajectories and their corresponding mean buoyancy across three sensitivity configurations: the Normal, WBF-off, and no-ice seeding cases (Figs. 9 and 12). To provide a more quantitative assessment of the time evolution, we analyze the accumulated altitude differences between these configurations in 100 s intervals (Fig. 11). During the initial stage of evolution (t<300 s), the trajectories remain nearly identical, which is physically explained by the evolution of plume buoyancy shown in Fig. 12. In the first 200 s, the buoyancy across all three ensembles shows minimal deviation, and the accumulated height differences remain close to zero (Fig. 11). This indicates that neither the initial latent heat of freezing nor the WBF growth exerts a significant immediate influence on the plume's altitude. During this period, vertical motion is almost entirely governed by the pre-existing environmental vertical velocity.

A distinct thermodynamic response emerges as the plumes evolve. Between t=200 and t=500 s, the buoyancy of both the Normal and WBF-off cases begins to increase significantly relative to the no-ice-nucleation baseline. As shown in Fig. 11, the height difference between the WBF-off and No-ice cases (green line) rises more rapidly than the difference between the Normal and WBF-off cases (blue line) during this interval. This suggests that the initial buoyancy gain is primarily attributed to the latent heat released during the freezing of cloud droplets and the onset of depositional growth at the condensation rate (L_c). This pulse of energy provides the impulse for the “secondary ascent” observed after 300 s.

The mature stage (t>600 s) is characterized by a secondary divergence in buoyancy, where the Normal plumes (solid lines in Fig. 12) exhibit a consistently higher rate of buoyancy increase compared to the WBF-off (L_s=L_c) plumes (dashed lines). Also, the accumulated height difference between the Normal and WBF-off cases (blue line in Fig. 11) shows a steady, non-linear increase over time. These two divergences represent the combined thermodynamic contribution of the WBF process and the enhanced depositional growth driven by adiabatic cooling. By extracting the latent heat of freezing (L_f) during the transition from both liquid and vapor phases to ice, these processes act as a persistent “internal engine” that sustains the plume's upward momentum.

This late-stage effect is quantified in Fig. 10 by calculating the ratio R_WBF:

$$\begin{array}{}\text{(2)}& R_{\text{WBF}}={\displaystyle \frac{Z_{\text{Normal}}-Z_{\text{WBF off}}}{Z_{\text{Normal}}-Z_{\text{no ice}}}}\end{array}$$

The results show that the excess latent heat accounts for approximately 25 % to 30 % of the total seeding-induced altitude gain. This stabilization in the averaged ratio is kinematically constrained by the shallow cloud layer in our numerical setup, as many plumes reach the cloud top after t>600 s. In individual cases where the plume remains within the cloud, the relative importance of the WBF process continues to increase over time. This ∼ 30 % contribution stems from two sources: (1) The complete elimination of the net latent heat gain from the WBF process (where droplet evaporation and ice growth now yield zero net energy change); (2) The reduction in latent heat released during depositional growth, as the vapor-to-ice transition is restricted to the lower condensation rate (L_c) instead of the sublimation rate (L_s).

Figure 10 The ratio of the mean-altitude difference between the Normal and WBF-off simulations to that between the Normal and No-ice seeding simulations at (a) 444 and (b) 800 s. Here, R_WBF is define in Eq. (2).

While the majority of the altitude gain (∼ 70 %) is driven by the initial freezing and the “baseline” latent heat of deposition (calculated at L_c), the divergence between the Normal and WBF-off plumes after t>600 s (Fig. 9) widens over time. This trend suggests a potential hierarchical partitioning in the thermodynamic drivers. We hypothesize that in the early stages of ice crystal growth, the buoyancy is primarily sustained by the rapid relaxation of ambient water vapor toward ice saturation – a process that occurs even in the reduced-heating case (L_s=L_c), despite lower energy efficiency. As the ice crystals grow larger and their total surface area increases, the relative importance of the WBF process (which relies on the liquid water reservoir) intensifies, acting as a secondary sustaining engine. Thus, while the ice-specific feedbacks appear secondary in the timeframe analyzed here, they may play a progressively dominant role in the long-term vertical transport and persistence of the seeded ice plume.

Figure 11 Time evolution of the accumulated plume-height differences between the three model configurations: Normal, WBF-off, and No-ice seeding. The height difference is calculated as the difference in mean plume altitude between two configurations within 100 s time intervals and averaged over all trajectories. Positive values indicate that the first configuration exhibits a higher plume altitude than the second configuration. Shaded regions represent the standard error of the mean across trajectories.

Figure 12 Same as Fig. 10, but for the time evolution of the mean buoyancy of the ice plume.

4.4 Microphysical evolution of the ice plume

The properties of the trajectories discussed in Sect. 4.2 influence the microphysical evolution of the seeded plume. Figure 13 presents the temporal evolution of the bulk ice properties. ICNC exhibits a rapid, monotonic decay across all simulations (Fig. 13a), driven primarily by turbulent dispersion and dilution rather than aggregation or sedimentation loss. While the decay rates are generally comparable, the downdraft-seeded plumes (blue lines) display a slightly more pronounced reduction in ICNC between t=100 and 400 s compared to updraft and no-wind cases. This aligns with the findings in Fig. 8, where downdraft plumes undergo greater vertical expansion, leading to stronger volumetric dilution. However, by t=800 s, ICNC values across all regimes converge, suggesting that dispersion eventually homogenizes the number density regardless of the initial perturbation.

In contrast to the convergence of number concentration, the evolution of ice particle size and mass shows significant divergence based on the vertical wind regime. The mean ice diameter (D_ice, Fig. 13b) increases continuously due to vapor deposition. Plumes in updrafts develop the largest crystals, followed by no-wind and downdraft cases. This stratification is even more distinct in the Ice Water Content (IWC, Fig. 13c). The updraft plumes, which are lofted rapidly into the supersaturated upper cloud regions, experience optimal growth conditions, resulting in the highest IWC. Conversely, downdraft plumes, which initially subside towards the cloud base or into sub-saturated air, experience suppressed growth rates, yielding significantly lower ice mass despite having similar number concentrations.

Figure 13 Temporal evolution of mean ice microphysical properties averaged within the plume. (a) Ice Crystal Number Concentration (ICNC), showing dilution due to dispersion. (b) Mean ice diameter (D_ice). (c) Ice Water Content (IWC). Shaded areas represent the standard deviation for all wind regimes.

Figure 14 Evolution of Ice Particle Size Distributions (PSD) sampled at downstream distances of (a) y=2, (b) y=3, and (c) y=4 km. Shaded areas are as in Fig. 13.

The evolution of the size distribution of the ice crystals confirms that while the growth rate varies, the growth mode remains consistent. Figure 14 displays the PSD at downstream distances of $y=\mathrm{2},\mathrm{3},$ and 4 km. Across all wind conditions, the PSDs exhibit a coherent rightward shift (towards larger diameters) over time. The PSD in the downdraft regime exhibits a broader (i.e., flatter) shape than in the other two regimes. This can be attributed to the larger vertical extent of the ice plume, which exposes ice particles to varying temperatures and, consequently, different diffusional growth rates at different altitudes. In addition, the initial downward trajectory exposes the ice particles to an environment with lower ice supersaturation, and part of the ice plume even descends below the cloud layer (see Appendix A). As a result, the number concentration of small particles (diameter <200 µm) remains relatively high, even after the ice plume has traveled 4 km.

The structural similarity of the PSDs in the no-wind and updraft regimes suggests that the seeding mechanism primarily determines the initial spectral characteristics, whereas the environmental updraft velocity mainly modulates the rate of spectral shift rather than altering the fundamental shape of the distribution. This is likely because, in both regimes, the ice plume follows a similar trajectory toward the cloud top, with the key difference being that the onset of upward motion occurs earlier in the updraft regime.

Figure 15 Same as Fig.  14, except that the PSD is shown at fixed altitudes of 350 (a–c), 450 (d–f), and 550 m (g–i).

Figure 16 Temporal response of the liquid phase within the ice plume. (a) Cloud Droplet Number Concentration (CDNC). (b) Mean liquid droplet diameter (D_liq). (c) Relative Humidity (RH) with respect to water. Shaded areas are as in Fig. 13.

The PSDs across the three regimes become increasingly similar at the same altitude in the simulations as the altitude increases (Fig. 15). This relative resilience of the PSD shape likely stems from the horizontally homogeneous quasi-uniform supersaturation environment within the liquid-dominated cloud layer, where ice crystals grow under similar vapor pressure deficits relative to water saturation regardless of their vertical trajectory. Furthermore, sedimentation does not yet play a significant role.

The growth of ice occurs at the expense of the liquid phase via the WBF process. This interaction is evident in Fig. 16, where CDNC decreases as the plume evolves. An anomaly is observed in the downdraft and no-wind simulations around t≈200 s: both CDNC (Fig. 16a) and mean liquid diameter (D_liq, Fig. 16b) exhibit a sharp minimum, coinciding with a drop in Relative Humidity (RH) (Fig. 16c). This “dip” corresponds to the period when subsiding plumes approach or penetrate cloud base, encountering drier air and causing partial evaporation or rapid glaciation. As the secondary lofting mechanism (driven by latent heat) takes over after t>300 s, the plumes re-ascend into the mixed-phase layer, allowing RH and liquid properties to stabilize. This highlights that seeding into downdrafts carries a risk of transient sublimation, although the plume ultimately recovers.

5 Conclusions and outlook

In this study, we investigated the dynamical and microphysical life cycle of glaciogenic ice plumes in supercooled stratiform clouds using the detailed bin microphysics model SCALE-AMPS. By constraining the simulations with in situ observations from the CLOUDLAB campaign, we successfully reproduced key plume characteristics – including ice number concentration, liquid water depletion, and particle size distributions – thereby establishing a robust framework for isolating the impact of environmental vertical wind conditions on seeding efficacy.

Our sensitivity analysis reveals that the evolution of a seeded ice plume is not merely a passive response to the initial advective field but is fundamentally shaped by latent heat release. The key findings are summarized as follows:

1. Thermodynamic decoupling and self-lofting: While the initial trajectory of the ice plume is kinematically governed by the vertical velocity at the injection moment, a decisive transition occurs after approximately 300–500 s. Regardless of whether the plume is seeded into an updraft or a downdraft, it eventually exhibits a mostly coherent net ascent. We attribute this universal lofting to the buoyancy generated by latent heat release from vapor deposition, which overcomes initial subsidence.

2. Evolution of the WBF Contribution: The relative importance of the WBF process – characterized by the net latent heat gain from the liquid-to-ice transition – increases over time. While the ensemble-averaged additional latent heat release specifically attributed to the WBF mechanism (L_s−L_c) stabilizes at approximately 25 %–30 % of the altitude gain after the first 15 min, this stabilization is primarily a consequence of the limited vertical extent of the simulated cloud layer, which kinematically constrains further plume ascent. In specific cases where the plume has not yet reached the cloud upper boundary, the relative contribution of the WBF mechanism continues to increase after t>600 s. We hypothesize a hierarchical growth preference: ice crystals first rapidly deplete the ambient vapor supersaturation (the primary early-stage driver) until the environment approaches ice saturation (e_i≈1), followed by a sustained growth phase increasingly supported by the evaporation of liquid droplets (the WBF process). This suggests that the WBF mechanism may become more important in the later stages of the plume's evolution as the available vapor reservoir is exhausted.

3. Vertical dispersion and detectability trade-off: The interaction between this induced ascent and the cloud-top inversion creates a trade-off. Plumes in updrafts ascend rapidly but are terminated by the inversion layer (“ceiling effect”). Conversely, plumes in downdrafts traverse the boundary layer more slowly, resulting in a larger vertical dispersion (up to ∼250 m spread). Crucially, this enhanced vertical dispersion compensates for the lower mean altitude of downdraft plumes, ensuring that the seeding signature remains geometrically detectable by downstream in situ sampling.

4. Microphysical resilience and transient desiccation: The size of the ice crystals increases consistently, and the size distribution is similar across all regimes. However, the persistence of the liquid phase is sensitive to the trajectory. Downdraft-seeded plumes exhibit a transient “dip” in liquid water content and relative humidity around t≈200 s as they approach the drier air below cloud base. Despite this temporary desiccation, the subsequent latent-heat-driven re-ascent allows the plume to recover and remain within the mixed-phase regime.

These findings have significant implications for both observational strategies and numerical modeling of precipitation enhancement. First, the identification of the “self-lofting” mechanism suggests that glaciogenic seeding is more dynamically robust than previously assumed; the success of seeding is not strictly contingent on targeting pre-existing updrafts, as the plume is capable of generating its own buoyancy to sustain its advection time in the mixed-phase layer. Second, the structural analysis confirms that Eulerian observations (e.g., fixed ground sites) are likely to capture the seeding signal, even in downdraft conditions, due to the compensatory effect of enhanced vertical dispersion.

Future work should expand this analysis to include a wider range of environmental conditions. Specifically, it would be interesting to vary the vertical wind shear profiles, as the strength of the shear directly modulates the boundary layer mixing processes and the lateral spreading of the ice plume. Additionally, sensitivity studies regarding thermodynamic profiles, such as temperature and latent water content (LWC), are warranted. These variables directly govern ice crystal habit regimes and riming efficiency, thereby influencing the magnitude of latent heat release.

Furthermore, the updrafts and downdrafts in this idealized framework are generated internally by boundary-layer mixing and therefore have characteristic horizontal scales of approximately 𝒪(100) m. In realistic atmospheric conditions, however, vertical motions in stratiform clouds also include larger-scale components, such as those associated with orographically forced gravity waves and mesoscale circulations. Incorporating vertical motions across a wider range of spatial scales would allow for a more comprehensive assessment of their interactions with smaller-scale motions and the associated microphysical – thermodynamic feedbacks.

Finally, since this study relied on a specific bin microphysics scheme, inter-comparison with bulk microphysics schemes would be valuable to determine whether standard parameterizations can adequately resolve the latent heat feedbacks that drive the secondary plume ascent. Ultimately, recognizing the active thermodynamic role of seeded ice plumes is essential for refining the representation of aerosol-cloud interactions for weather modification simulations.