Secondary ice production affects tropical convective clouds under different aerosol conditions

Secondary ice production (SIP) modulates tropical convection, but its interactions with aerosol loading remain uncertain. Using the UK Met Office Unified Model with detailed SIP parameterizations, this study investigates how multiple SIP mechanisms affect a Hector storm under weak instability and varying cloud condensation nuclei (CCN) concentrations. A reference simulation including all SIP mechanisms at a prescribed droplet number concentration (N_d) of 400 cm⁻³ is first evaluated against satellite, radar and aircraft observations. At moderate aerosol loading (N_d = 400 cm⁻³), SIP substantially enhances the ice phase. Ice number increases by over 3 orders of magnitude at temperatures warmer than −10 °C, and ice water content rises by up to ∼ 120 % in the anvil. SIP also reduces outgoing longwave radiation (OLR) by 6.8 W m⁻² (−2.5 %), increases outgoing shortwave radiation (OSR) by 7.9 W m⁻² (+2 %), over a 110 km × 110 km region in 6 h. For this case, precipitation becomes more organized near the convective core, with total accumulation increasing by ∼ 25 %. Under polluted (N_d = 800 cm⁻³) and clean (N_d = 200 cm⁻³) conditions, SIP effects are weaker or spatially fragmented, indicating non-linear aerosol dependence. Sensitivity experiments show that, when combined with rime-splintering (RS), Mode 1 (droplet fragmentation during freezing) weakens at high CCN owing to reduced supercooled raindrop supply in the mixed-phase layer, while Mode 2 (drop–ice collisions with splashing) is inhibited at low CCN due to fewer encounters with large rimed ice particles. In contrast, RS and its combination with ice–ice collisional breakup show relatively limited aerosol sensitivity in this case.

1 Introduction

Tropical deep convection plays a central role in the global climate and hydrological cycle by regulating precipitation, radiation, and energy transport. It also redistributes aerosols and trace gases, shapes the vertical structure of temperature and humidity in the troposphere, and contributes to setting the tropical tropopause through radiative–convective equilibrium (e.g., Houze, 2004; Fueglistaler et al., 2009; Sherwood et al., 2014). The microphysical processes governing the development of these systems remain a major source of uncertainty in weather and climate models (e.g., Bony et al., 2015; Grabowski, 2015; Morrison et al., 2020). In particular, the formation of ice in mixed-phase clouds is not fully explained by primary ice nucleation alone, as observed ice crystal number concentrations (ICNC) often exceed measured ice-nucleating particle (INP) concentrations by several orders of magnitude (e.g., Hobbs and Rangno, 1985; Ladino et al., 2017; Field et al., 2017). This discrepancy points to the importance of secondary ice production (SIP), whereby new ice particles are generated from pre-existing ones.

A range of SIP mechanisms has been identified, including rime-splintering (hereafter RS; Hallett and Mossop, 1974), collisional breakup of ice particles (hereafter BR; Vardiman, 1978; Grzegorczyk et al., 2023), and sublimation-induced fragmentation (Oraltay and Hallett, 1989). Droplet shattering during freezing has also been documented and is commonly classified into Mode 1 and Mode 2 based on the symmetry of freezing (Dye and Hobbs, 1968; Lauber et al., 2018; Phillips et al., 2018; James et al., 2021). Laboratory experiments confirm that these processes can generate large numbers of new ice particles, while SIP signatures have also been reported in both mid-latitude and tropical convection, as reviewed by Field et al. (2017), Korolev and Leisner (2020) and Huang et al. (2022). At the microphysical level, SIP substantially enhances ICNC and ice water content (IWC), and alters the phase distribution between liquid and ice (Phillips et al., 2017a; James et al., 2023; Han et al., 2024).

At the storm scale, modeling and observational studies have shown that SIP can strongly influence cloud radiative and precipitation responses. By increasing anvil extent and optical depth, SIP modifies outgoing longwave radiation (OLR) and shortwave reflection (Young et al., 2019; Hawker et al., 2021a). At the same time, changes in riming and melting pathways modify rainfall intensity and organization. For instance, Han et al. (2024) reported that including SIP in central European convection increased ICNC by more than twentyfold and reduced surface precipitation by up to 20 %. Grzegorczyk et al. (2025a) found that enabling SIP in a heavy-precipitation case improved agreement with in-situ ice measurements while reducing heavy rain occurrence. Dedekind et al. (2021) showed that in mid-latitude winter convection, SIP enhanced ICNC and surface precipitation, with impacts highly sensitive to the chosen parameterizations. In tropical convection, SIP has been associated with the maintenance of extensive anvils and changes in precipitation efficiency (Hawker et al., 2021b; Qu et al., 2022; Waman et al., 2022), while in stratiform and Arctic mixed-phase clouds it helps explain persistent ICNC well above INP limits (Sotiropoulou et al., 2020; Zhao et al., 2021). Despite these advances, substantial uncertainties remain. Different SIP mechanisms dominate in different temperature ranges, their relative importance varies across environments, and parameterizations diverge among models (Field et al., 2017; Sullivan et al., 2018; Korolev and Leisner, 2020).

Beyond SIP itself, changes in CCN concentration are expected to modulate its efficiency by reshaping hydrometeor properties in the mixed-phase region. Under lower CCN, collision–coalescence tends to occur earlier, accelerating warm-rain processes. This can increase the abundance of raindrops that participate in riming, but it also depletes the supercooled liquid aloft that sustains SIP (Rosenfeld and Woodley, 2000; Khain et al., 2005). With increasing cloud condensation nuclei (CCN), more numerous, smaller droplets delay autoconversion and loft more liquid above the freezing level, maintaining supercooled water and modifying droplet shattering and collisional breakup between ice particles (Rosenfeld et al., 2008; Qu et al., 2020). Which SIP route dominates therefore depends on temperature, the intensity of riming, and the availability of precipitation hydrometeors suitable for riming, implying a non-linear dependence of SIP efficiency on aerosol loading (Fan et al., 2016; Grabowski and Morrison, 2020; Varble et al., 2023). Qu et al. (2020) showed that aerosol-controlled droplet spectra modulate SIP and the partition between liquid and ice, influencing storm evolution. Previous studies have also linked CCN impacts to changes in SIP and subsequent storm electrification (Phillips and Patade, 2022; Huang et al., 2025). However, most existing studies examine aerosol effects and SIP sensitivities separately. Systematic evaluations of aerosol–SIP interactions across CCN regimes, particularly in tropical convection, remain limited (e.g., Field et al., 2017; Korolev and Leisner, 2020).

This study addresses this gap by investigating how SIP interacts with CCN to influence cloud–precipitation–radiation processes in tropical convection. We focus on a weak-CAPE Darwin storm on 6 February 2006 and conduct sensitivity experiments with the Unified Model and CASIM microphysics across clean, moderate, and polluted CCN levels. We first evaluate a reference simulation of the complete model, including all SIP mechanisms, against the available observations, and then use matched sensitivity experiments to diagnose the effects of SIP and CCN relative to this reference state. In addition, we examine the CCN sensitivity of RS and RS-based configurations with additional SIP mechanisms (Mode 1, Mode 2, and BR), to clarify their relative contributions. Through this analysis, we aim to provide new insights into how SIP and aerosols jointly modulate tropical convection and storm evolution, and to emphasize the importance of representing SIP processes accurately in weather and climate models.

2 Method

2.1 Model description

The Met Office Unified Model (UM) vn13.0 with the Cloud AeroSol Interacting Microphysics Module (CASIM; Field et al., 2023) is used in this study. Simulations are conducted using the Regional Atmosphere and Land (RAL) 3.2 configuration (Bush et al., 2023). The UM employs the ENDGame non-hydrostatic dynamical core with a semi-Lagrangian formulation (Wood et al., 2014). Subgrid cloud fraction is diagnosed using the bimodal cloud fraction scheme of Van Weverberg et al. (2021) for coupling to radiation in the RAL3.2 configuration, and no convection parameterization is active. Subgrid mixing in the boundary layer is represented by the Lock et al. (2000) scheme. Further details are provided in Bush et al. (2025).

Cloud microphysics are represented using the CASIM scheme, configured in two-moment mode to prognostically predict both mass and number concentrations for five hydrometeor species: cloud droplets, raindrops, ice crystals, graupel, and snow. Particle size distributions follow a generalized gamma distribution. In CASIM, the cloud droplet number concentration (N_d) can be either prescribed as a fixed value or diagnosed from background aerosol. To systematically examine aerosol sensitivity and isolate the effects of prescribed droplet-number conditions, we adopt the fixed N_d approach in this study. As in RAL3.2, a vertical taper is imposed on the fixed N_d to avoid unrealistically high droplet concentrations at upper levels and the associated excessive homogeneous freezing (Bush et al., 2025).

Rain is generated through auto-conversion, accretion, and melting of ice-phase particles, while losses include evaporation and freezing by ice. The separate treatment of mass and number concentrations allows for variability in hydrometeor size and fall speed, which directly affects sedimentation and evaporation. In contrast to single-moment schemes where particle size is fixed for a given mass, CASIM allows dynamic adjustment of size distributions based on aerosol and microphysical evolution.

Heterogeneous ice nucleation is parameterized using the default scheme of Cooper (1986), which prescribes ice crystal number concentrations as a function of temperature. Immersion freezing of raindrops is parameterized following Bigg (1953), with freezing probability depending on drop volume and temperature. When freezing occurs, rain mass and number are converted to graupel, providing a bulk representation of raindrop freezing in CASIM (Field et al., 2023; Miltenberger et al., 2020). In addition, rain–ice interactions are represented through accretion/collection processes in CASIM, and rain collecting ice can contribute to snow or graupel formation depending on the hydrometeor state and model configuration. Homogeneous freezing of cloud droplets is parameterized to occur below −38 °C (Field et al., 2023). In this study, aerosol perturbations are represented by prescribed N_d. INP concentrations are not varied and are not linked to N_d; therefore coupled CCN–INP effects are not represented. In CASIM, aggregation and accretion among ice-phase hydrometeors are represented through parameterizations for self-collection and inter-species collection processes involving ice, snow, and graupel. The associated collection efficiencies are prescribed in the model and are temperature-dependent.

In addition to primary ice nucleation, this study incorporates multiple secondary ice production (SIP) mechanisms. The default SIP process in CASIM is rime-splintering (RS), which produces ice splinters through the riming of supercooled droplets by graupel and snow between −2.5 and −7.5 °C, with peak efficiency at −5 °C (Hallett and Mossop, 1974; Mossop, 1985). The production rate is based on the cloud liquid accreted by graupel and snow, using a default efficiency of 350 ice splinters per milligram of rimed liquid. In our previous study (Sun et al., 2025), we implemented three additional SIP mechanisms into CASIM: ice-ice collisional breakup (BR), and droplet shattering during symmetrical and asymmetrical freezing (Mode 1 and Mode 2). The BR parameterization follows the theoretical formulation of Phillips et al. (2017b), and accounts for collisions between ice-phase particles (graupel, snow, and ice crystals). The production rate depends on collision kinetic energy, rimed fraction, temperature, and particle size. For droplet shattering, we adopt the two-mode parameterization from Phillips et al. (2018). Mode 1 (M1) describes the collisions between frozen drops and smaller ice particles, with splinter numbers fitted to laboratory data using Lorentzian curves, peaking around −15 °C. Mode 2 (M2) involves collisions between freezing raindrops and more massive ice particles, where fragmentation is governed by the ratio of collision kinetic to surface energy. An empirical factor of 0.3, representing the ice-containing fraction of fragments, is applied following James et al. (2021).

Radiative transfer is handled by the SOCRATES scheme (Edwards and Slingo, 1996; Manners et al., 2024), which employs the two-stream approximation for fluxes and spherical harmonics for radiance. The single-scattering properties of ice crystals are derived from both the ice water content and the number concentration of ice-phase hydrometeors predicted by CASIM. This allows changes in ice particle characteristics associated with different SIP mechanisms to influence the simulated cloud optical properties. Turbulence is represented by the blended boundary layer scheme (Boutle et al., 2014), which incorporates both local and non-local mixing processes.

2.2 Model set-up and design of the simulations

Simulations were conducted for the 6 February 2006 case from the Aerosol and Chemical Transport in tropIcal convection (ACTIVE) campaign (Vaughan et al., 2008; Connolly et al., 2013). The ACTIVE campaign was carried out in the Darwin region of northern Australia and aimed to investigate aerosol-cloud interactions in tropical convection. This case, which occurred during the break period of the monsoon season, featured a relatively weaker convective storm over the Tiwi Islands (Fig. 1a). In contrast to the vigorous pre-monsoon convection, the environment was characterized by a dry mid-troposphere and limited convective instability, with a convective available potential energy (CAPE) of approximately 1010 J kg⁻¹ (Fig. 1b). The sounding also suggests a warm cloud base and a relatively deep warm-cloud layer, which likely favor warm-rain formation during the early stage of storm development.

Figure 1 (a) Overview of model orography (terrain height, unit: m) and spatial domains used for radar and satellite data analysis (adapted from Sun et al., 2025). The full simulation domain is a 1350 km × 1350 km square centered at (12.4° S, 130.8° E). The subdomain used for radar evaluation and model diagnostics is indicated by the white box. The effective range of the CPOL radar, centered at (12.25° S, 131.05° E), is shown as the overlaid circle. (b) Atmospheric profiles from Darwin Airport station (12.42° S, 130.89° E) at 00:00 UTC on 6 February 2006. The green and red solid lines represent dew point and temperature, and the blue line represents the parcel lapse rate.

The cloud-resolving simulations are conducted at a horizontal resolution of 1.5 km, using a domain of 900 × 900 grid points centered over the Darwin region (12.4° S, 130.8° E). The vertical grid contains 90 levels, extending up to 40 km, with finer resolution in the lower atmosphere (52 levels below 10 km, 33 below 4 km, and 16 below 1 km). Each simulation starts at 18:00 UTC on 5 February 2006 and runs for 36 h. Boundary and initial conditions are prescribed from the ECMWF Reanalysis version 5 (ERA5) reanalysis data, updated every 6 h. A spin-up period of 6 h is applied, with a model time step of 75 s. In this configuration, the semi-Lagrangian dynamical core permits relatively large time steps. To avoid hydrometeors crossing multiple vertical levels within a single time step, sedimentation is treated using an approach based on Rotstayn (1997), as implemented in the UM (Wilson and Ballard, 1999; Field et al., 2023). To better represent surface features, we incorporate the Shuttle Radar Topography Mission (SRTM; NASA, 2013) dataset as ancillary input. For analysis, we focus on a 1° × 1° subdomain over the Tiwi Islands (indicated by the white box in Fig. 1a), selected to minimize potential boundary effects.

Table 1 Model experiment overview.

Note: RS refers to rime-splintering (Hallett and Mossop, 1974); M1 and M2 indicate two droplet shattering modes (Phillips et al., 2018); BR denotes ice-ice collisional breakup (Phillips et al., 2017b).

To investigate the sensitivity of mixed-phase convection to secondary ice production (SIP) mechanisms and aerosol conditions, a series of simulations was conducted by combining three prescribed in-cloud droplet number concentrations (N_d) with different SIP configurations (Table 1). A fixed N_d approach was employed to provide controlled droplet-number conditions and to allow the microphysical, radiative, and precipitation responses to SIP to be compared consistently across three different aerosol environments. The rates of the different SIP pathways are linked to N_d through the prescribed droplet field and its influence on warm-rain production, providing the basis for assessing SIP impacts across the three prescribed N_d conditions. N_d was specified at 200, 400, and 800 cm⁻³, corresponding to low, intermediate, and high droplet-number conditions (hereafter referred to as “clean”, “moderate”, and “polluted”, respectively). These values span the range reported during ACTIVE, where cloud droplet number concentrations of approximately 150, 500, and 900 cm⁻³ were associated with the monsoon-break, pre-monsoon, and burning-period environments, respectively (Connolly et al., 2013). The three prescribed values therefore provide a consistent framework for examining how the SIP response varies across plausible aerosol conditions.

For each N_d value, six SIP scenarios were tested: all SIP mechanisms active (allSIP), no SIP process enabled (noSIP), rime-splintering only (RS), and three combined schemes involving RS with one additional mechanism: RS with Mode 1 (RS+M1), RS with Mode 2 (RS+M2), and RS with ice–ice collisional breakup (RS+BR). For each prescribed N_d, the noSIP simulation at the same N_d provides the matched comparison for the corresponding SIP experiments, so that SIP effects are diagnosed with N_d held fixed. The matched comparison between allSIP and noSIP is used to quantify the combined effect of SIP, while the RS, RS+M1, RS+M2, and RS+BR simulations are used to examine the contributions and aerosol dependence of individual mechanisms. RS was retained as the baseline SIP configuration because it is the default SIP process in CASIM. This design enables a systematic assessment of how aerosol conditions influence the impacts of SIP on cloud microphysics, radiation, and precipitation. Of these simulations, allSIP-400 is used as the reference simulation, that is, the complete model configuration, with all SIP mechanisms active, that provides the most observationally defensible representation of the observed storm. Its prescribed N_d is selected by comparing the three complete-model configurations (allSIP-200, allSIP-400 and allSIP-800) with the satellite, radar and aircraft observations (Sect. 3.1). Although the ACTIVE monsoon-break value (∼ 150 cm⁻³; Connolly et al., 2013) lies closer to the clean configuration, allSIP-400 agrees more closely with the observations than allSIP-200, and the uncertainty in N_d is represented by the 200 and 800 cm⁻³ experiments. The remaining simulations are treated as perturbations of the SIP treatment and of N_d (Sect. 3.2–3.4).

2.3 Observational data

The modeled reflectivity fields were evaluated against observations from the C-band polarimetric radar (CPOL), located at Gunn Point (12.25° S, 131.05° E). The CPOL radar provides volumetric scans every 10 min, with a vertical resolution of 0.5 km and horizontal resolution of 2.5 km. In the model, radar reflectivity was diagnosed from the simulated hydrometeor particle size distributions using a forward operator based on a simple Rayleigh scattering assumption, consistent with Field et al. (2023). To facilitate comparison with the model, we used the Statistical Coverage Product (May et al., 2008), which represents the time–height distribution of the fraction of a defined domain exceeding 10 dBZ, a threshold chosen to capture the spatial and temporal characteristics of convective precipitation. The same 10 dBZ threshold was applied to the simulated reflectivity fields for comparison. The analysis domain is indicated in Fig. 1a.

Simulated rainfall accumulations were also compared against CPOL radar–derived precipitation datasets, available at 10 min intervals and 1 km horizontal resolution. Although residual uncertainty may arise from reflectivity calibration and C-band attenuation, CPOL quantitative precipitation estimates generally provide a reliable representation of tropical rainfall (Louf et al., 2019; Jackson et al., 2021). Aircraft observations from the ACTIVE campaign were also considered, with the NERC Dornier sampling the lower troposphere up to approximately 4 km and the Egrett sampling the upper-level anvil outflow at approximately 12–14 km (Vaughan et al., 2008). Available in situ measurements of IWC and ice particle number concentration were confined to upper-level clouds and did not directly sample the mixed-phase region where SIP is expected to be most active. These upper-level observations are compared with the reference simulation in Fig. 4. Radiosonde profiles from Darwin Airport station (12.42° S, 130.89° E) were used to characterize the thermodynamic environment. Additionally, top-of-atmosphere radiation was evaluated using equivalent outgoing longwave radiation (OLR) calculated from Multifunctional Transport Satellite (MTSAT) brightness temperature, with spatial and temporal resolutions of 1.25 km and 30 min, respectively. The observed and simulated OLR are compared in Fig. 2.

Figure 2 Evaluation of the reference simulation (allSIP-400) against MTSAT-derived equivalent outgoing longwave radiation (OLR; W m⁻²). (a–c) Observed OLR at 07:30, 08:30 and 09:30 UTC on 6 February 2006. (d–f) Simulated OLR from allSIP-400 at 06:30, 07:30 and 08:30 UTC; simulation times are shifted by −1 h to match the earlier convective development in the model. (g) Domain-averaged OLR from allSIP-400 (red line) and MTSAT (grey markers; error bars show ±1 standard error).

3 Results

Section 3.1 first evaluates the reference simulation (allSIP-400) against the available observations. Section 3.2–3.4 then analyse the sensitivity experiments, in which SIP effects are diagnosed from matched SIP–noSIP differences at fixed N_d, and the dependence on N_d is assessed across the three prescribed values. Observations are not used in these sensitivity sections.

3.1 Evaluation of the reference simulation

On 6 February 2006, the Hector storm developed over the Tiwi Islands during the break period of the monsoon season. MTSAT imagery (Fig. 2a–c) shows shallow cumulus initiating along the sea-breeze front (Connolly et al., 2013), developing rapidly between 07:30 and 08:30 UTC, and reaching maturity by 09:30 UTC. The coldest cloud tops, with OLR below 150 W m⁻², were concentrated over the central island, with an expanding westward anvil. With the simulation time shifted by −1 h to account for the earlier development of the simulated storm, allSIP-400 reproduces the location of the storm and its westward anvil expansion (Fig. 2d–f), although its minimum OLR remains higher than observed. In the domain mean (Fig. 2g), allSIP-400 captures the observed decrease and recovery of OLR. Its mature-stage minimum is lower than observed by ∼ 7 W m⁻², and OLR is underestimated before and after the mature stage by up to ∼ 15 W m⁻².

Figure 3 Evaluation of the reference simulation (allSIP-400) against CPOL radar observations. Time–height CFADs of reflectivity greater than 10 dBZ from (a) CPOL and (b) allSIP-400. Accumulated precipitation (mm) during 04:00–10:00 UTC on 6 February 2006 from (c) CPOL and (d) allSIP-400. Domain-averaged (e) precipitation rate (mm h⁻¹) and (f) accumulated precipitation (mm) from CPOL (grey) and allSIP-400 (red); the grey dashed line in panel (e) is a fitted curve to the observations. Simulation times are shifted by −1 h. The 0, −20, −40 and −60 °C isotherms are shown as dashed lines in panel (b).

The CPOL reflectivity CFAD (Fig. 3a) shows echoes exceeding 10 dBZ extending to about 15 km during the mature stage. allSIP-400 reproduces this vertical extent, with frequent echoes above 10 km (Fig. 3b). The simulated rainband structure and heavy rainfall cores (> 100 mm) are broadly consistent with the CPOL accumulation (Fig. 3c, d). The peak domain-mean precipitation rate in allSIP-400 (∼ 1.35 mm h⁻¹) is close in timing and magnitude to the observed peak (Fig. 3e). Accumulated precipitation is underestimated, at ∼ 3.54 mm compared with ∼ 4.8 mm observed (Fig. 3f).

Figure 4 Vertical profiles of (a) ice number concentration (N_ice) and (b) ice water content (IWC) from the reference simulation (allSIP-400; red) and the corresponding simulation without SIP (noSIP-400; blue), compared with Egrett in situ observations. Shading shows the model interquartile range. Both are calculated over regions where IWC > 0.01 g m⁻³ during 02:30–08:30 UTC on 6 February 2006. Observations are grouped into 0.5 km height bins: black dots show bin means, horizontal bars the first and third quartiles, and numbers the sample count in each bin. Grey dashed lines indicate the 0, −20, −40 and −60 °C isotherms. The aircraft sampled only the upper-level anvil and not the mixed-phase region.

Egrett observations sampled the anvil outflow at about 12–14 km. Within this layer, the ice particle number concentration and IWC of allSIP-400 are broadly consistent with the observations (Fig. 4; noSIP-400 is also shown, to indicate how far the SIP response extends into the sampled layer). These observations do not sample the mixed-phase region between about 5 and 9 km, where SIP is most active. However, the upper-level ice population reflects the cumulative history of the convective ascent, including the supercooled liquid lofted to homogeneous-freezing levels, which is itself influenced by SIP lower in the cloud (Phillips and Patade, 2022). The other two complete-model configurations agree less well with these observations, which is the basis for choosing allSIP-400 as the reference. allSIP-200 produces a higher minimum OLR (∼ 165 W m⁻²), a lower peak precipitation rate (∼ 0.95 mm h⁻¹) and a smaller accumulation (∼ 3.25 mm). allSIP-800 produces a lower peak precipitation rate (∼ 0.94 mm h⁻¹), a smaller accumulation (∼ 3.11 mm), and a more fragmented anvil and rainfall pattern than observed.

Overall, the reference simulation reproduces the main observed features of the storm: the timing and location of convection, the depth of the radar echoes, the organisation and peak intensity of rainfall, and the upper-level ice properties. The main remaining differences are the OLR minimum, the OLR before and after the mature stage, the accumulated rainfall, and the absence of in situ observations in the mixed-phase layer. These differences provide the context for interpreting the sensitivity experiments that follow.

3.2 SIP impacts on radiation under varying CCN conditions

Figure 5 shows simulated OLR from the allSIP simulations at the three prescribed N_d values (allSIP-200, allSIP-400 and allSIP-800), with the corresponding noSIP fields in Fig. S1 in the Supplement. The underlying microphysical mechanisms are discussed in Sect. 3.4.1.

Figure 5 Spatial distribution of simulated outgoing longwave radiation (OLR; W m⁻²) from allSIP simulations for (a–c) N_d = 200 cm⁻³, (d–f) N_d = 400 cm⁻³ and (g–i) N_d = 800 cm⁻³ at 06:30, 07:30 and 08:30 UTC, respectively. N_d denotes the prescribed in-cloud droplet number concentration. Simulation times are shifted by −1 h to match the earlier convection development in the model. Equivalent plots for simulated noSIP cases are provided in Fig. S1 in the Supplement.

The impacts of SIP on radiation were first examined by comparing the spatial distribution of OLR between each allSIP simulation and its corresponding noSIP case at 07:30 UTC (Fig. S1). For N_d = 200 cm⁻³, allSIP-200 shows a weaker OLR reduction relative to noSIP-200 (Figs. 5b, S1b), with a minimum OLR of ∼ 165 W m⁻² and a smaller increase in the low-OLR coverage (∼ 10 %–15 %). For N_d = 400 cm⁻³, allSIP-400 produces a minimum OLR of ∼ 150 W m⁻² (Fig. 5e), about 25 W m⁻² lower than noSIP-400 (∼ 175 W m⁻²), and expands the area with OLR < 200 W m⁻² by around 25 %–30 % (Fig. S1e). At N_d = 800 cm⁻³, allSIP-800 still exhibits reduced OLR relative to noSIP-800, but the low-OLR region becomes more fragmented (Figs. 5h, S1h).

Figure 6 Time-height plots of CFADs (contoured frequency by altitude diagrams) for radar reflectivity greater than 10 dBZ from (a) allSIP-200, (b) noSIP-200, (c) allSIP-400, (d) noSIP-400, (e) allSIP-800 and (f) noSIP-800 simulations. The 0, −20, −40, and −60 °C isotherms are indicated by dashed lines in all panels.

The CFADs (contoured frequency by altitude diagrams) of reflectivity (> 10 dBZ), as shown in Fig. 6, further link these radiative patterns to storm structure. For N_d = 200 cm⁻³, the difference between allSIP-200 and noSIP-200 is smaller and mainly confined to mid-levels (Fig. 6a, b), consistent with the weaker radiative impact (Fig. 5a–c). For the N_d = 400 cm⁻³ pair, allSIP-400 sustains a deeper vertical extent of reflectivity > 10 dBZ above 10 km than noSIP-400 (Fig. 6c, d), consistent with the expanded low-OLR region. At N_d = 800 cm⁻³, reflectivity > 10 dBZ reaches upper levels (Fig. 6e, f), but the echoes are more vertically dispersed and occur with lower frequency at each height than in the N_d = 400 cm⁻³ pair. This suggests a deeper yet optically thinner distribution of hydrometeors, which may explain the less concentrated low-OLR region (Fig. 5g–i). Overall, the intermediate N_d simulations produce the most persistent upper-level reflectivity, particularly when all SIP mechanisms are included. The response to SIP is weaker under low N_d, while the high N_d simulations exhibit a more vertically dispersed reflectivity structure.

Figure 7 Time series of domain-averaged (a) outgoing longwave radiation (OLR) and (b) outgoing shortwave radiation (OSR) from the simulations. Blue and red lines denote noSIP and allSIP simulations, respectively. Dotted lines denote N_d = 200 cm⁻³, solid lines denote N_d = 400 cm⁻³, and solid lines with square markers denote N_d = 800 cm⁻³. The same time adjustment is applied to all simulations to account for the earlier development of the simulated storm.

Table 2 Comparison of radiation responses to SIP across CCN levels.

Note: All values are calculated over the analysis region shown in Fig. 1a, during 02:30–08:30 UTC on 6 February 2006. ΔOLR and ΔOSR represent the differences between SIP and noSIP simulations. Relative change (%) is calculated as (OLR_SIP − OLR_noSIP) $/$ OLR_noSIP × 100 %, and similarly for OSR.

The temporal evolution of domain-averaged OLR is shown in Fig. 7a, and the SIP-induced changes are summarized in Table 2. In all simulations, OLR decreases after approximately 02:00 UTC as convection develops, reaches a minimum during the mature stage, and recovers towards the end of the analysis period. At N_d = 400 cm⁻³, allSIP-400 reaches a substantially lower minimum than noSIP-400. The SIP-induced change is smaller at N_d = 800 cm⁻³, and small and of opposite sign at N_d = 200 cm⁻³. Quantitatively, the strongest reduction occurs at N_d = 400 cm⁻³, where the allSIP case drops below 250 W m⁻² and the mean allSIP–noSIP difference reaches −6.86 W m⁻² (−2.5 %). At N_d = 200 cm⁻³, allSIP-200 produces a slight increase of +0.61 W m⁻² (+0.2 %) relative to noSIP-200, whereas at N_d = 800 cm⁻³, allSIP-800 produces a decrease of −3.34 W m⁻² (−1.2 %) relative to noSIP-800. The reductions under the intermediate and high N_d conditions are consistent with stronger anvil development and enhanced upper-level ice content.

In contrast, outgoing shortwave radiation (OSR; Fig. 7b) increases after 02:00 UTC and peaks around 06:00 UTC, but its response to SIP is more variable. At N_d = 400 cm⁻³, SIP enhances OSR relative to noSIP, with peak values exceeding 500 W m⁻² compared to ∼ 480 W m⁻², and a mean increase of +7.95 W m⁻² (+2.0 %). At N_d = 200 cm⁻³, however, SIP reduces OSR by −6.55 W m⁻² (−1.8 %), while at N_d = 800 cm⁻³ the reduction is smaller at −4.05 W m⁻² (−0.9 %). This contrast suggests that SIP has a robust and consistent impact on longwave radiation by maintaining colder, more extensive anvils, whereas its shortwave effect is weaker and dependent on CCN concentration. In addition, OSR differences may be further influenced by variations in cloud phase, particle size, and the timing of solar insolation (Yin and Porporato, 2017).

3.3 SIP impacts on precipitation under varying CCN conditions

The spatial distribution of accumulated precipitation from 04:00–10:00 UTC (Fig. 8) shows that SIP affects both the intensity and organization of rainfall, with responses depending on CCN concentrations. At N_d = 200 cm⁻³ (Fig. 8a, b), allSIP-200 slightly increases the number of rainfall cores relative to noSIP-200 but does not substantially change the overall coverage, indicating a weaker response in this case. allSIP-400 produces a more organized rainband with more heavy rainfall cores (> 100 mm) than noSIP-400 (Fig. 8c, d), indicating stronger organization when SIP is included. Under N_d = 800 cm⁻³ (Fig. 8e–f), allSIP-800 broadens the precipitation area relative to noSIP-800, but the rainfall pattern is more fragmented than in allSIP-400, resembling the more dispersed pattern in the OLR distribution (Fig. 5h–i).

Figure 8 Spatial distribution of accumulated precipitation (mm) during 04:00–10:00 UTC on 6 February 2006 from (a) allSIP-200, (b) noSIP-200, (c) allSIP-400, (d) noSIP-400, (e) allSIP-800 and (f) noSIP-800.

Figure 9 Time series of domain-averaged (a) precipitation rate (mm h⁻¹) and (b) accumulated precipitation (mm) from simulations without SIP (noSIP; blue) and with all SIP mechanisms activated (allSIP; red) for N_d = 200 cm⁻³ (dotted lines), N_d = 400 cm⁻³ (solid lines), and N_d = 800 cm⁻³ (solid lines with square markers). The time axis of the simulations has been shifted by −1 h, as in Fig. 3.

Domain-averaged precipitation time series (Fig. 9) further quantify these effects. The peak precipitation rate in allSIP-400 (∼ 1.35 mm h⁻¹) is notably higher than in noSIP-400 (∼ 0.80 mm h⁻¹; Fig. 9a). In the N_d = 400 cm⁻³ pair, SIP enhances the production of ice-phase particles aloft, accelerating mixed-phase processes (Sun et al., 2025). This promotes melting and sustains a robust rain region, as shown in Fig. 8c. For N_d = 200 cm⁻³, allSIP-200 (∼ 0.95 mm h⁻¹) produces a lower peak rate than noSIP-200 (∼ 1.20 mm h⁻¹). In clean conditions, early precipitation formation likely involves a stronger contribution from warm-rain processes. SIP accelerates glaciation during the early storm stage, shifting part of the condensate from liquid to ice phase and likely reducing the liquid water available for warm-rain collision–coalescence. This is consistent with the lower peak precipitation rates under clean conditions. The precipitation duration remains comparable, resulting in similar accumulated totals between allSIP-200 (∼ 3.25 mm) and noSIP-200 (∼ 3.27 mm; Fig. 9b). For N_d = 400 cm⁻³, allSIP-400 (∼ 3.54 mm) yields a higher accumulated total than noSIP-400 (∼ 2.81 mm, +25 %; Fig. 9b), consistent with its stronger precipitation peak. At N_d = 800 cm⁻³, high aerosol loading delays precipitation onset due to smaller cloud droplet sizes and reduced collision–coalescence efficiency. In this case, SIP increases ice water content, leading to greater late-stage precipitation in allSIP-800 than in noSIP-800 (∼ 2.30 mm). However, the total accumulation in allSIP-800 (∼ 3.11 mm) remains lower than in allSIP-400 (∼ 3.54 mm), as the limited warm-rain formation under high droplet number concentrations cannot be fully compensated by the enhanced ice-phase process. The simulated domain totals span approximately 50 % across the combinations of SIP and N_d (2.30–3.54 mm), demonstrating the sensitivity of precipitation accumulation to both prescribed droplet number and SIP configuration. This provides a clear motivation to represent and constrain these processes more accurately in precipitation prediction.

Overall, these results suggest that SIP can intensify rainfall and produce a more concentrated precipitation pattern, with the strongest response in allSIP-400 relative to noSIP-400. At N_d = 200 cm⁻³, large cloud droplets allow warm-rain processes to occur efficiently even without SIP, so the additional ice-phase pathway has only a small impact on precipitation. At N_d = 800 cm⁻³, enhanced cloud droplet numbers delay warm-rain processes and lead to more dispersed precipitation; allSIP-800 increases totals relative to noSIP-800 but does not reach the allSIP-400 total.

3.4 SIP impacts on cloud microphysics under varying CCN conditions

3.4.1 Microphysics associated with radiation and precipitation responses

To interpret the radiative and precipitation responses presented in previous sections, we examine the vertical structures of key hydrometeors under different CCN and SIP conditions (Fig. 10). These include cloud droplets, raindrops, ice crystals, graupel, total ice water content (IWC), and updraft strength. Cloud droplet number concentration increases with CCN as expected, ranging from ∼ 0.5 × 10⁸ m⁻³ at N_d = 200 cm⁻³ to over 2.0 × 10⁸ m⁻³ at N_d = 800 cm⁻³ (Fig. 10a), while the mean-mass radius of cloud droplets decreases with increasing CCN (Fig. 10f). The impact of SIP on cloud droplets is negligible compared to the dominant effect of CCN.

Figure 10 Vertical profiles of horizontally averaged microphysical and dynamical properties from simulations without SIP (noSIP) and with all SIP mechanisms (allSIP), under three CCN concentrations (N_d = 200, 400, and 800 cm⁻³). Panels (a)–(d) show number concentrations of cloud droplets, rain, ice crystals, and graupel, respectively. Panels (f)–(i) show the corresponding mean-mass radii. Panels (e) and (j) display the ice water content (IWC, sum of ice crystals, graupel, and snow) and 95th percentile updraft velocity (W_max), used as a proxy of maximum updraft strength. Panels (a)–(b) and (f)–(g) are averaged over regions where either the ice water path or rain water path exceeds 1 g m⁻², while panels (c)–(e) and (h)–(j) are averaged over regions where IWC > 0.01 g m⁻³ during 02:30–08:30 UTC on 6 February 2006. In panels (a)–(e) and (j), dotted lines denote N_d = 200 cm⁻³, solid lines denote N_d = 400 cm⁻³, and solid lines with square markers denote N_d = 800 cm⁻³. Gray dashed lines indicate the 0, −20, −40, and −60 °C isotherms.

SIP alters the vertical distribution of raindrops, and its impact varies with CCN levels. Raindrop number concentration (N_rain) decreases from the 200 to 800 cm⁻³experiments (Fig. 10b), consistent with reduced warm-rain efficiency at higher CCN. Near the melting level (∼ 4–5 km), N_rain is controlled primarily by CCN and differs little between allSIP and noSIP at a given CCN. Below 4 km, however, N_rain is generally higher in allSIP-200, allSIP-400, and allSIP-800 than in their noSIP counterparts, indicating that melting of SIP-generated ice contributes additional raindrops to lower levels. Above 5 km, SIP increases N_rain at N_d = 200 cm⁻³ but decreases it at N_d = 400 and 800 cm⁻³. In the clean pair, allSIP-200 shows higher N_rain than noSIP-200 at 5–7 km, suggesting earlier melting of the relatively limited ice produced aloft. By contrast, allSIP-400 and allSIP-800 both show lower N_rain than noSIP-400 and noSIP-800 above 5 km, implying that more condensate is retained in frozen form aloft before being converted to rain at lower levels. This behaviour suggests a stronger ice-phase contribution to rainfall in allSIP-400 and allSIP-800 than in the clean pair, as further examined using the ice and IWC diagnostics below.

For raindrop size, the mean-mass radius below 5 km is similar across CCN levels and between noSIP and allSIP runs (Fig. 10g), indicating that mixing, coalescence, and evaporation tend to homogenize raindrop sizes at low levels. Above 5 km, raindrops are smallest at N_d = 200 cm⁻³ and larger at N_d = 400 and 800 cm⁻³. In particular, allSIP-200 shows a smaller radius than noSIP-200 at 5–7 km, with allSIP-200 reaching < 160 µm compared with ∼ 200 µm in noSIP-200. This likely reflects a higher N_rain with little change in rainwater content, reducing the mean mass per drop and thus the radius. In contrast, allSIP-400 and allSIP-800 produce larger raindrops than noSIP-400 and noSIP-800 by ∼ 30–100 µm above 5 km, consistent with a contribution from melting of enhanced ice-phase hydrometeors to fewer but larger drops.

Time–height sections of rainwater content (RWC) further support this pathway (Fig. S2 in the Supplement). Compared to noSIP-400, allSIP-400 shows a pronounced RWC enhancement at ∼ 3–6 km from ∼ 04:00–06:00 UTC (Fig. S2b), followed by a near-surface (0–2 km) RWC peak around ∼ 06:30 UTC that coincides with the precipitation-rate maximum in Fig. 9a. The enhancement is weaker and delayed in allSIP-800 relative to allSIP-400, consistent with suppressed warm-rain development at high CCN. Nevertheless, allSIP-800 still exceeds noSIP-800 (Fig. S2c, f, i), indicating that SIP continues to enhance rainwater under polluted conditions through the melting of SIP-enhanced ice-phase hydrometeors. In the clean pair, differences between allSIP-200 and noSIP-200 are small or negative, especially between 04:00–06:00 UTC (Fig. S2g), indicating that SIP only weakly modifies rainwater mass under low-CCN conditions. Taken together, these rain diagnostics are consistent with the contrasting precipitation responses across the three CCN regimes: the clearest enhancement occurs in allSIP-400, the response in allSIP-800 is weaker and delayed, and differences between allSIP-200 and noSIP-200 remain small.

SIP leads to substantial enhancement in ice crystal number concentration (N_ice) between 5 and 9 km, especially in allSIP-400 (Fig. 10c). At 5–7 km, allSIP-400 reaches ∼ 3 × 10³ m⁻³, more than three orders of magnitude higher than noSIP-400. AllSIP-800 shows a similar enhancement in this layer, but the increase is more vertically dispersed, with a secondary peak (∼ 0.5 × 10³ m⁻³) near 8 km. In the clean case, the N_ice enhancement in allSIP-200 is confined below 7 km, and values between 7–9 km (∼ 0.2 × 10³ m⁻³) are roughly half those in allSIP-400 and allSIP-800. The mean-mass radius of ice crystals (Fig. 10h) also increases in allSIP runs, especially in allSIP-400, where values exceed 100 µm near 7 km compared with ∼ 60 µm in noSIP-400. In contrast, the noSIP cases show little variation in radius across CCN levels between 5–9 km. Previous studies have shown that SIP tends to produce numerous smaller ice particles with reduced sedimentation (e.g. Dedekind et al., 2021). In this case, however, the enhanced crystal population is associated with larger mean-mass radii in allSIP-400 and, to a lesser extent, in allSIP-800. This suggests that part of the SIP-generated crystal population continues to grow after formation, likely through vapor deposition aided by stronger updrafts and greater IWC aloft (Field et al., 2023). In the clean pair, the radius difference between allSIP-200 and noSIP-200 remains small, likely because both updraft and IWC changes are small under low-CCN conditions.

Graupel number concentration (N_graupel) is also enhanced by SIP across all CCN levels (Fig. 10d). The largest increases occur in allSIP-200 and allSIP-400, both of which reach ∼ 10³ m⁻³ near 6 km, compared to ∼ 10 m⁻³ in the corresponding noSIP runs. The enhancement is smaller in allSIP-800, consistent with previous findings that polluted conditions suppress graupel formation by delaying raindrop freezing (Sun et al., 2021, 2024). In contrast to the ice crystals, graupel radius decreases in all SIP runs (Fig. 10i). For example, at ∼ 5.5 km, graupel radius is ∼ 300 µm in allSIP-200 and ∼ 400 µm in noSIP-200. This suggests that SIP increases the number of graupel particles, while the available supercooled liquid is distributed among more particles, limiting the growth of individual graupel through riming. By contrast, the diagnosed snow category shows a weaker response, suggesting that the SIP signal in snow may be partly moderated by aggregation and sedimentation in this case (Field et al., 2023). This is also consistent with a previous Hector simulation with CASIM (e.g., Sun et al., 2025). We therefore focus on crystals and graupel here.

Table 3 Comparison of IWC and updraft responses to SIP mechanisms across CCN levels.

Note: ΔIWC is taken at the time of maximum difference between SIP and noSIP simulations. Relative change (%) is calculated as (IWC_SIP − IWC_noSIP) $/$ IWC_noSIP × 100 %. Peak altitudes correspond to the vertical levels where ΔIWC or W_max reaches its maximum for each configuration.

Together, increases in ice crystal and graupel enhance total ice water content (IWC) and the associated latent heating. The strongest IWC enhancement occurs in allSIP-400, where values exceed 0.35 g m⁻³ above 7 km, compared with ∼ 0.25 g m⁻³ in noSIP-400 at the same altitude (Fig. 10e). Near 10 km, allSIP-400 shows a maximum enhancement of +121.7 % relative to noSIP (Table 3). The response in allSIP-800 is weaker and more vertically dispersed, with peaks near 7 and 10 km and a maximum difference of +0.16 g m⁻³ (+87.9 %) around 10 km. In the clean pair, the IWC difference between allSIP-200 and noSIP-200 remains small despite increases in both N_ice and N_graupel. Under low-CCN conditions, larger cloud droplets and earlier rainfall suggest a stronger warm-rain contribution, which likely reduces available supercooled liquid content. At the same time, the relatively dry environment and weak updrafts are less favorable for sustained vapor deposition, limiting further ice-mass growth. As a result, the IWC difference remains small, and the associated allSIP–noSIP updraft differences are also weak (Fig. 10j). By contrast, allSIP-400 and allSIP-800 produce larger and more persistent IWC aloft (Fig. 10e). The latent heating from freezing strengthens mid- to upper-level updrafts by 2–4 m s⁻¹ relative to noSIP-400 and noSIP-800 (Fig. 10j), which helps maintain condensate aloft and supports further ice growth.

Figure 11 Time–height cross sections of ice water content (IWC; g m⁻³) for simulations with all SIP mechanisms included: (a) allSIP-200, (b) allSIP-400, and (c) allSIP-800; and without SIP: (d) noSIP-200, (e) noSIP-400, and (f) noSIP-800. Panels (g)–(i) show the differences between allSIP and noSIP simulations at each CCN concentration (N_d = 200, 400, 800 cm⁻³; i.e., allSIP minus noSIP). Panels (a)–(f) are averaged over regions where IWC > 0.01 g m⁻³. The 0, −20, −40, and −60 °C isotherms are shown by the dashed lines. Corresponding plots for rain water content (RWC) are shown in Fig. S2 in the Supplement.

Time–height cross sections of IWC (Fig. 11) further illustrate how SIP modifies vertical cloud structure across aerosol conditions. IWC increases occur primarily above the melting level (∼ 5 km) during 04:00–08:30 UTC (Fig. 11a–f). The response is strongest at N_d = 400 cm⁻³, with deeper and longer-lasting IWC enhancement at upper levels than in the corresponding noSIP run (Fig. 11b, e, h), consistent with the largest IWC increase in Fig. 10e. At N_d = 800 cm⁻³, enhancement is weaker and delayed relative to N_d = 400 cm⁻³ (Fig. 11c, f, i), which is consistent with the weaker OLR suppression. At N_d = 200 cm⁻³, the changes induced by SIP are small and short-lived (Fig. 11a, d, g), indicating limited mass growth in the clean case. The difference panels (Fig. 11g–i) show that the SIP-induced IWC increase is concentrated at mid- to upper levels (∼ 7–13 km), with the strongest enhancement in the 400 cm⁻³ case (Fig. 11h). Together with the rain and RWC diagnostics discussed above, these differences help explain why the radiative response is strongest and the precipitation enhancement is clearest in allSIP-400.

3.4.2 Sensitivity to individual SIP mechanisms

We examine how different secondary ice production pathways, namely allSIP, RS, RSM1, RSM2, and RSBR, modify the storm relative to the corresponding noSIP simulations at N_d = 200, 400, and 800 cm⁻³, as shown in Figs. 10 and 12. Diagnostics are expressed as height-dependent enhancements over noSIP (e.g., ΔIWC = IWC_SIP − IWC_noSIP; see Table 3).

At the case level with all mechanisms active, SIP substantially increases IWC, with the largest enhancement at N_d = 400 cm⁻³. The maximum ΔIWC reaches about +0.19 g m⁻³ (+121.7 %), and the IWC maximum is close to 10 km (Fig. 10e and Table 3). The peak updraft shows the same nonmonotonic dependence on CCN. W_max attains 12.6 m s⁻¹ at N_d = 400 cm⁻³, compared with 9.1 and 10.5 m s⁻¹ at N_d = 200 and 800 cm⁻³ respectively (Fig. 10j). The altitude of W_max is higher at N_d = 400 cm⁻³, around 5.5 to 7 km, and lower in the clean and polluted cases (Table 3).

Figure 12 Vertical profiles of horizontally averaged (a–d) ice water content (IWC) and (e–h) 95th percentile updraft velocity (W_max) from simulations without SIP (noSIP) and with different SIP configurations: (a, e) RS, (b, f) RSM1, (c, g) RSM2, and (d, h) RSBR. Dotted lines denote N_d = 200 cm⁻³, solid lines denote N_d = 400 cm⁻³, and solid lines with square markers denote N_d = 800 cm⁻³. All profiles are averaged over regions with IWC > 0.01 g m⁻³ during 02:30–08:30 UTC on 6 February 2006. Gray dashed lines indicate the 0, −20, −40, and −60 °C isotherms.

The CCN sensitivity of ΔIWC depends on the mechanism. RSM1 peaks at N_d = 400 cm⁻³ with +0.12 g m⁻³ (+78.1 %), remains moderate at N_d = 200 cm⁻³ with +0.10 g m⁻³ (+47.4 %), and becomes weak at N_d = 800 cm⁻³ with +0.03 g m⁻³ (+14.1 %; Fig. 12b). RSM2 also peaks at N_d = 400 cm⁻³ with +0.17 g m⁻³ (+109.7 %), reduces to +0.08 g m⁻³ (+44.1 %) at N_d = 800 cm⁻³, and is weakest at N_d = 200 cm⁻³ with +0.03 g m⁻³ (+16.4 %; Fig. 12c). RS and RSBR produce large enhancements at N_d = 200 and 400 cm⁻³, approximately +0.13–0.16 g m⁻³ (+85 %–106 %), but weaker responses at N_d = 800 cm⁻³, approximately +0.06–0.08 g m⁻³ (+29 %–40 %; Fig. 12a, d, and Table 3).

CCN also modulates the altitude of the SIP-induced maxima (Fig. 12 and Table 3). For RSM1, the altitude of the ΔIWC maximum shifts from about 10 km at N_d ≤ 400 cm⁻³ to about 6 km at N_d = 800 cm⁻³. For RSM2, the maximum remains near 10 km at N_d ≥ 400 cm⁻³ but lowers to about 7 km at N_d = 200 cm⁻³. For RS and RSBR, the peak rises from about 7 km at N_d = 200 cm⁻³ to about 10 km at N_d = 400 cm⁻³. These vertical shifts are consistent with the stronger upper-level ice in the N_d = 400 cm⁻³ case and the weaker enhancement under polluted conditions, which helps explain the CCN dependence of the OLR reduction. Peak updraft response is weaker than the IWC response, with ΔW_max rarely exceeding 1.5 m s⁻¹, and the largest W_max occurring at RS-400 and RS-800, RSM1-200, RSM2-800, and RSBR-400 (Fig. 12e–h).

Figure 13 Vertical profiles of horizontally averaged (a–d) ice crystal number concentration (N_ice) and (e–h) ice crystal mean-mass radius from simulations without SIP (noSIP) and with different SIP configurations: (a, e) RS, (b, f) RSM1, (c, g) RSM2, and (d, h) RSBR. Dotted lines denote N_d = 200 cm⁻³, solid lines denote N_d = 400 cm⁻³, and solid lines with square markers denote N_d = 800 cm⁻³. All profiles are averaged over regions with IWC > 0.01 g m⁻³ during 02:30–08:30 UTC on 6 February 2006. Gray dashed lines indicate the 0, −20, −40, and −60 °C isotherms. Equivalent plots for cloud droplet and rain properties are shown in Figs. S3–S4 in the Supplement.

Microphysical diagnostics provide further detail. As shown in Fig. 13, for RS and RSBR, ice crystal number concentrations exhibit evident peaks (> 10³ m⁻³) between 5 and 8 km, and the three CCN profiles are nearly overlapping, indicating limited aerosol sensitivity. In the same layer, ice crystal mean-mass radii range from ∼ 50–120 µm and show little variation with CCN (Fig. 13e, h). For RSM1, the three CCN profiles are close together between 5 and 7 km, while between 7 and 9 km the N_d = 800 cm⁻³ case is clearly weaker than the N_d = 200 and 400 cm⁻³ cases. RSM1 produces the largest crystal radii among the mechanisms, often ∼ 150–200 µm between 5 and 9 km, but these radii are nearly invariant with CCN. To test whether the suppression of RSM1 under polluted conditions arises from limited raindrop supply, we examine vertical profiles of cloud droplet and raindrop properties (Figs. S3–S4 in the Supplement), focusing on number concentrations and mean-mass radii from about 0 to −20 °C, corresponding to about 5–9 km in these simulations. As shown in Fig. S3, cloud droplets are primarily controlled by CCN rather than by the SIP configuration. Number concentrations increase from N_d = 200 to 800 cm⁻³, whereas mean-mass radii decrease with CCN. In contrast, differences between RSM1 and noSIP are minor compared with the CCN effect. The prevalence of numerous but smaller droplets under high CCN (Fig. S3) is consistent with slower collision–coalescence and delayed growth toward raindrop sizes, implying a source limitation for Mode 1 in the 0 to −20 °C layer. Raindrop fields show the corresponding signature for RSM1 (Fig. S4). Between 5–7 km, raindrop number concentration decreases with CCN, with the ordering N_d = 200 > 400 > 800 cm⁻³, while the radius is largest at N_d = 400 cm⁻³, smaller at N_d = 800 cm⁻³, and smallest at N_d = 200 cm⁻³. Between 7–9 km, the number concentration retains the same ordering, and the radius at N_d = 400 and 800 cm⁻³ is similar and both exceed that at N_d = 200 cm⁻³ (Fig. S4). The larger radii above 5 km at N_d = 400 and 800 cm⁻³ are consistent with the concentration of available liquid into fewer drops at higher levels and may be aided by the melting of ice-phase hydrometeors. Taken together, these patterns indicate that polluted conditions do not favor the availability of the supercooled raindrops required by Mode 1. Although raindrops are similar or larger at higher levels for N_d = 800 cm⁻³, their number and residence time in this layer are reduced, so the frequency of Mode 1 is lower. The reduced raindrop supply at N_d = 800 cm⁻³ also aligns with the weaker RSM1 ice production between 7–9 km (Fig. 13) and with the minimum ΔIWC reported for RSM1 at high CCN (Table 3). Collectively, these diagnostics support the interpretation that the suppression of RSM1 under polluted conditions arises from a warm-rain source limitation.

Figure 14 Same as Fig. 13, but for (a–d) graupel number concentration (N_graupel), and (e–h) graupel mean-mass radius.

For RSM2, the ice crystal number concentration shows little sensitivity to CCN below ∼ 7 km, whereas between 7 and 9 km it is highest at N_d = 400 cm⁻³, lower at N_d = 800 cm⁻³, and lowest at N_d = 200 cm⁻³ (Fig. 13c). Unlike Mode 1 (RSM1), Mode 2 (RSM2) relies on collisions between raindrops and large ice particles, such as graupel. Therefore, we examine graupel characteristics to assess the efficiency of Mode 2 under different CCN levels. As shown in Fig. 14, graupel number in 5–9 km exceeds noSIP and generally decreases with increasing CCN. The graupel radius is smallest at N_d = 200 cm⁻³ (∼ 240 µm) and larger at 400–800 cm⁻³ (∼ 400 µm) between 5 and 7 km, and this contrast remains at higher levels. Clean conditions (N_d = 200 cm⁻³) also promote earlier warm rain, reducing supercooled liquid at upper levels and limiting liquid–ice overlap with sufficiently large graupel. Together, smaller graupel and reduced overlap lower the likelihood of effective collisions and suppress Mode 2 efficiency. In contrast, at N_d = 400 cm⁻³, both liquid water and graupel size are more favorable, consistent with the ΔIWC maximum at N_d = 400 cm⁻³ and the minimum at N_d = 200 cm⁻³. Overall, Mode 2 is suppressed under clean conditions due to smaller graupel, weaker liquid–ice overlap, and the reduced likelihood of effective collisions.

4 Discussion

4.1 SIP–CCN effects on radiation and precipitation

The experiments show that SIP alters the radiative and precipitation characteristics of the 6 February 2006 storm, but its impacts depend on CCN concentrations. Under moderate aerosol loading (N_d = 400 cm⁻³), SIP produces the largest reduction in outgoing longwave radiation (OLR), expansion of low-OLR anvil areas, and the strongest enhancement of organized precipitation. By contrast, SIP effects are weaker under clean (N_d = 200 cm⁻³) and polluted conditions (N_d = 800 cm⁻³), and the OSR response is overall modest and varies in sign with CCN. These results highlight a non-linear SIP–CCN influence on storm structure and radiative forcing.

Agreement between the reference simulation and the available observations (Sect. 3.1) does not exclude compensating errors. The observations constrain the storm-scale radiative, reflectivity and precipitation evolution and the upper-level ice field, but not the mixed-phase layer in which SIP operates. Offsetting biases in, for example, primary ice nucleation, raindrop freezing or aggregation therefore cannot be ruled out. The SIP responses diagnosed from the matched experiments should thus be interpreted as model responses relative to an observationally consistent reference state, rather than as validated estimates of SIP in the observed storm.

The underlying mechanisms can be linked to the modulation of cloud microphysics by CCN. At low CCN, large droplets form and precipitation develops early through warm-rain processes, depleting supercooled liquid water and leaving little available for SIP-induced ice growth. As a result, OLR reduction and precipitation changes are limited. At high CCN, numerous small droplets delay warm-rain processes and sustain supercooled liquid water at higher altitudes. The reduced collision–coalescence efficiency leads to the lowest accumulated rainfall in noSIP-800 among the three CCN levels (Fig. 9b). When SIP is activated, this additional supercooled liquid water provides a source for secondary ice production. In this case, the extra ice particles enhance riming and graupel formation, and subsequent melting increases surface rainfall, so allSIP-800 yields more rainfall than noSIP-800. The SIP-related increase remains smaller than at N_d = 400 cm⁻³, because at moderate CCN the mixed-phase layer contains a more favorable combination of supercooled liquid water, raindrops, and large rimed particles within the temperature ranges where the SIP mechanisms are active, allowing more efficient ice growth and melt out. At moderate CCN, the SIP mechanisms therefore act together most effectively, leading to stronger upper-level ice growth, thicker and more organized anvils, a larger reduction in OLR (Figs. 5, 7), and enhanced rainfall (Figs. 8, 9). Previous studies also suggest that SIP impacts on precipitation can depend strongly on the parameterization used, and may increase ice particle number while suppressing localized or heavy precipitation (e.g., Dedekind et al., 2021; Grzegorczyk et al., 2025b). Thus, SIP partly compensates for warm rain suppression under polluted conditions, but its strongest combined radiative and precipitation impact occurs at moderate CCN, where the microphysical environment most strongly supports secondary ice production and conversion to surface rainfall.

This non-linear response should also be interpreted in the context of the weak-CAPE environment of this case (CAPE ≈ 1010 J kg⁻¹; Fig. 1b). With limited instability, the balance between warm-rain and ice-phase processes could be more sensitive to CCN perturbations. In contrast, our earlier pre-monsoon case (Sun et al., 2025), characterized by higher CAPE and stronger convection, showed larger reductions in OLR, stronger increases in OSR, and suppressed precipitation in allSIP compared to noSIP. This enhanced radiative signal arose from deeper and more persistent anvils formed under high CAPE. Vigorous updrafts transported large amounts of hydrometeors into the upper troposphere, strengthening both longwave cooling and shortwave reflection (Feng et al., 2011; Hartmann and Berry, 2017). At the same time, stronger glaciation converted more supercooled liquid water into ice particles, enhancing the formation of precipitation particles in the upper levels. However, much of this ice was retained in anvils or sublimates/melts before reaching the ground, leading to reduced surface rainfall despite greater ice production (Han et al., 2024; Grzegorczyk et al., 2025a). By contrast, in the weak-CAPE 6 February 2006 storm, OLR responses remain robust but OSR is only weakly affected (Figs. 5, 7). SIP enhances OSR at N_d = 400 cm⁻³, but this effect reverses at N_d = 200 cm⁻³ and nearly vanishes at N_d = 800 cm⁻³. This weak and non-monotonic OSR response likely reflects the combined influences of changes in anvil coverage, optical depth and particle size, and the timing of SIP-induced anvils relative to the maximum solar insolation. As a result, shortwave fluxes exhibit the SIP signal much less clearly than longwave cooling in this case. This behavior is consistent with Finney et al. (2025), who showed across 12 convective cases that changes in cloud droplet number, ice nucleating particles, and secondary ice configurations can modify anvil albedo and high cloud area, but that these shortwave responses are generally small compared with the case to case spread, highlighting the strong role of environmental variability and the need for caution when interpreting single case results. Meanwhile, precipitation increases in allSIP relative to noSIP (Figs. 8, 9). Under weak instability, cloud depth and anvil extent are limited compared to the high-CAPE case. SIP generates additional ice particles that likely promote riming and subsequent melting, which may provide an extra source of surface rainfall. This may help explain why SIP increases precipitation in the 6 February 2006 case, whereas the pre-monsoon case showed the opposite response.

These results underscore the importance of considering environmental conditions, such as CAPE, when assessing SIP-CCN interactions in tropical convection. Across both radiation and precipitation responses, the strongest enhancement consistently occurs under moderate CCN, indicating that SIP plays an important role in amplifying aerosol effects on tropical convection, particularly under weakly unstable environments.

4.2 Sensitivity to SIP mechanisms and model uncertainties

The aerosol response of secondary ice production is mechanism dependent and non-linear. RSM1 is suppressed under polluted conditions (N_d = 800 cm⁻³), RSM2 is inhibited under clean conditions (N_d = 200 cm⁻³), and RS together with RSBR show limited CCN sensitivity. The combined effect produces the strongest signals at moderate CCN (N_d = 400 cm⁻³), including enhanced upper-level ice, higher peak updrafts, and the largest reduction of outgoing longwave radiation (Figs. 12–14; Table 3).

For RSM1, the diagnostics indicate a limitation in the warm-branch supply. Mode 1 depends on the freezing of supercooled raindrops and their collisions with small ice particles. Under high CCN, cloud droplets are more numerous and smaller, collision–coalescence proceeds more slowly, and the delivery of raindrops into the 0 to −20 °C layer is reduced. Consistently, RSM1 exhibits less ice number between 7 and 9 km and the smallest ΔIWC at N_d = 800 cm⁻³, while crystal sizes remain nearly unchanged across CCN levels (Fig. 13; Table 3).

For RSM2, Mode 2 depends on encounters between raindrops and large rimed ice particles. In clean conditions, warm rain forms and precipitation begins earlier, which limits the overlap of supercooled liquid and large ice in the 7–9 km layer. The diagnostics show the lowest upper-level ice number and the smallest graupel radii at N_d = 200 cm⁻³, whereas moderate to high CCN maintains larger graupel and stronger upper-level ice production, with ΔIWC maximized at N_d = 400 cm⁻³ (Fig. 14; Table 3).

Our simulations further demonstrate that the allSIP configuration produces the strongest enhancement at moderate CCN (N_d = 400 cm⁻³). This non-linear behavior reflects the need for supercooled raindrops and sufficiently large rimed ice to co-exist within the mixed-phase layer. Under clean and polluted conditions, one of the droplet-shattering pathways is limited: Mode 1 at high CCN due to the reduced supply of supercooled raindrops, and Mode 2 at low CCN because collisions with large rimed ice particles are less frequent in the mixed-phase layer. At moderate CCN, supercooled liquid is available, and raindrops occur with larger rimed ice, allowing Mode 1 and Mode 2 to operate concurrently on top of the relatively CCN-insensitive RS and BR processes. As a result, the combined SIP contribution is largest at intermediate aerosol loading, yielding the greatest ΔIWC and clear dynamical and radiative responses (Figs. 7, 9, 10). This peak at intermediate aerosol loading is consistent with previous observational and modeling studies (Rosenfeld et al., 2008; Fan et al., 2016), and our results identify SIP as a process pathway that helps produce this non-linear CCN dependence. This response likely arises because the amounts of supercooled liquid water and rimed ice, which are controlled by N_d, do not change in a simple way across the temperature ranges where each SIP mechanism is most active. Different SIP parameterizations, with other temperature ranges or production rates, could change this relationship and in some cases produce a more monotonic response to N_d. However, the general result that aerosol control of liquid and ice in the mixed-phase layer affects the strength of SIP is expected to remain.

Beyond this CCN dependence, our results also indicate that droplet shattering mechanisms (Mode 1 and Mode 2) contribute more strongly to ice production than ice–ice collisional breakup (BR). This suggests that interactions between raindrops and ice particles can play an important role in sustaining secondary ice production under the simulated conditions. However, Grzegorczyk et al. (2025c) found that BR and RS were the most effective SIP mechanisms in reproducing ice numbers in deep convection, while droplet shattering played a minor role. Previous tagging-tracer analyses in continental deep convection also found a dominant contribution from BR to total ice number, whereas droplet shattering contributed comparatively less (e.g., Waman et al., 2022). Part of this difference may arise from differences in the microphysics schemes, as well as from the environmental conditions in our simulations, where graupel sizes and collision energies are relatively modest (Sun et al., 2025), limiting the efficiency of BR. One specific example is the treatment of raindrop freezing in the two modeling systems. UM-CASIM employs a formulation based on Bigg (1953) with modified parameter values, which is represented separately from the heterogeneous ice nucleation scheme and is not explicitly linked to INP activity. This treatment affects the number of freezing raindrops and graupel particles formed in the mixed phase region where droplet shattering can occur. This differs from recent approaches in which raindrop freezing is related to INP concentrations, such as Waman et al. (2022). In comparison, the bin microphysics model DESCAM uses the original Bigg (1953) scheme (Grzegorczyk et al., 2025c), which tends to produce fewer freezing raindrops (e.g., Morrison et al., 2005; Fan et al., 2009; Sullivan et al., 2018) and thus fewer opportunities for droplet shattering, so that BR and RS account for a larger share of the simulated secondary ice. A further uncertainty concerns the representation of heterogeneous ice nucleation. In the present study, this follows Cooper (1986) and may affect the relative enhancement attributed to SIP. The magnitude of the simulated SIP contribution should therefore be interpreted with caution. Environmental differences, such as CAPE and liquid water content, may also contribute to the SIP responses. These results highlight the uncertainties in quantifying the relative importance of SIP mechanisms. Consistent parameterizations and further observational constraints are still needed to separate the contributions of droplet shattering and collisional breakup in deep convective clouds.

5 Conclusions

This study examines how secondary ice production (SIP) interacts with aerosol loading to affect tropical convection, using the 6 February 2006 Darwin storm as a weak-CAPE case. The results show that SIP significantly alters cloud radiative properties, precipitation characteristics, and storm microphysics, with the strongest effects occurring at moderate aerosol concentrations (N_d = 400 cm⁻³). A reference simulation including all SIP mechanisms at N_d = 400 cm⁻³ reproduces the main observed features of the storm, including the timing and location of convection, the depth of the radar echoes, the peak rainfall rate and the upper-level ice properties, although it underestimates accumulated rainfall and does not reproduce the observed OLR evolution exactly.

Microphysically, SIP substantially enhances the ice phase. At N_d = 400 cm⁻³, ice crystal number concentration (N_ice) exceeds 3 × 10³ m⁻³ at temperatures warmer than −10 °C, more than three orders of magnitude above noSIP. Ice crystals also grow larger, with mean-mass radii exceeding 100 µm, and ice water content (IWC) increases by up to 121.7 %. The diagnosed snow response to SIP remains weak in this case, even though these changes could potentially favor conversion to snow. Relative to the corresponding noSIP runs, graupel number concentration increases at all CCN levels, although the graupel radius decreases due to enhanced competition for liquid water. These changes strengthen latent heating and vertical motion, with the strongest updrafts (W_max = 12.6 m s⁻¹) at N_d = 400 cm⁻³, compared with 10.5 m s⁻¹ at N_d = 800 cm⁻³ and 9.1 m s⁻¹ at N_d = 200 cm⁻³.

At the macroscale, SIP enhances both radiation and precipitation responses, but primarily under moderate aerosol loading. At N_d = 400 cm⁻³, SIP reduces outgoing longwave radiation (OLR) by 6.86 W m⁻² (−2.5 %) and expands the low-OLR cloud cover by ∼ 25 %–30 %. In contrast, the OLR response is weaker and spatially fragmented at N_d = 800 cm⁻³ (−3.34 W m⁻², −1.2 %), and small and of opposite sign at N_d = 200 cm⁻³ (+0.61 W m⁻², +0.2 %). Changes in outgoing shortwave radiation (OSR) are smaller and vary in sign, with a clear increase only at N_d = 400 cm⁻³ (+7.95 W m⁻²).

SIP also reorganizes rainfall patterns. At N_d = 400 cm⁻³, rainfall becomes more organized, producing stronger and more concentrated precipitation patterns. The domain-averaged peak rainfall rate rises from 0.80 to 1.35 mm h⁻¹, and total accumulation increases by ∼ 25 %. At N_d = 800 cm⁻³, rainfall becomes more spatially fragmented, and at N_d = 200 cm⁻³ the impact is small or negative, with peak rainfall rate suppressed.

Sensitivity experiments show distinct CCN responses among SIP mechanisms. RSM1 (RS + Mode 1) is most active at moderate CCN (ΔIWC = +0.12 g m⁻³ at N_d = 400 cm⁻³) but weakens at high CCN (+0.03 g m⁻³ at N_d = 800 cm⁻³), owing to a reduced supply of supercooled raindrops in the mixed-phase layer. RSM2 (RS + Mode 2) is strongest at moderate CCN (ΔIWC = +0.17 g m⁻³ at N_d = 400 cm⁻³), but weakens at lower CCN (+0.03 g m⁻³ at N_d = 200 cm⁻³), where collisions with larger rimed ice particles are less frequent. RS and RSBR are relatively less sensitive to CCN. Their combined contribution peaks at intermediate CCN, producing the largest ΔIWC and clear dynamical and radiative responses.

Nevertheless, the CCN sensitivity of individual SIP mechanisms remains uncertain, owing to differences in the adopted parameterizations, microphysical assumptions in the model, and environmental conditions. In addition, the physical understanding and representations of SIP processes are still subject to substantial uncertainty and ongoing debate. The available aircraft observations provide a useful constraint on upper anvil ice but do not directly constrain ice multiplication within the mixed-phase region. Future work should combine targeted modelling with co-located in situ aircraft measurements in the mixed-phase region, particularly of ice particle number concentration and particle size distributions, to evaluate the modelled ice-phase properties and constrain SIP responses across aerosol regimes.