Secondary ice production in summer clouds over the Antarctic coast: an underappreciated process in atmospheric models

The correct representation of Antarctic clouds in atmospheric models is crucial for accurate projections of the future Antarctic climate. This is particularly true for summer clouds which play a critical role in the surface melting of the ice shelves in the vicinity of the Weddell Sea. The pristine atmosphere over the Antarctic coast is characterized by low concentrations of ice nucleating particles (INPs) which often result in the formation of supercooled liquid clouds. However, when ice formation occurs, the ice crystal number concentrations (ICNCs) are substantially higher than those predicted by existing primary ice nucleation parameterizations. The rime-splintering mechanism, thought to be the dominant secondary ice production (SIP) mechanism at temperatures between −8 and −3 ^∘C, is also weak in the Weather and Research Forecasting model. Including a parameterization for SIP due to breakup (BR) from collisions between ice particles improves the ICNC representation in the modeled mixed-phase clouds, suggesting that BR could account for the enhanced ICNCs often found in Antarctic clouds. The model results indicate that a minimum concentration of about ∼ 0.1 L⁻¹ of primary ice crystals is necessary and sufficient to initiate significant breakup to explain the observations, while our findings show little sensitivity to increasing INPs. The BR mechanism is currently not represented in most weather prediction and climate models; including this process can have a significant impact on the Antarctic radiation budget.

1 Introduction

Predictions of Antarctic climate are hampered by the poor representation of mixed-phase clouds over the Southern Ocean and the Antarctic seas (Haynes et al., 2011; Flato et al., 2013; Bodas-Salcedo et al., 2014; Hyder et al., 2018). Model simulations reveal significant discrepancies in the Antarctic surface radiation budget associated with cloud biases that are driven by errors in the representation of the cloud microphysical structure (Lawson and Gettelman, 2014; King et al., 2015; Listowski and Lachlan-Cope, 2017). A correct representation of the cloud radiative impacts largely depends on the parameterization of cloud microphysical processes (Listowski and Lachlan-Cope, 2017; Hines et al., 2019; Young et al., 2019) and precipitation (Vignon et al., 2019) which determine the concentration and characteristics of liquid drops and ice crystals.

Ice crystals form at temperatures above −38 ^∘C through heterogeneous nucleation (Pruppacher and Klett, 1997); this means that the presence of insoluble aerosols that act as ice nucleating particles (INPs) is required. However, Antarctica and the Southern Ocean are relatively clean regions, and INPs are sparse (McCluskey et al., 2018; Schmale et al., 2019; Welti et al., 2020). Thus it is especially surprising that enhanced ice crystal number concentrations (ICNCs) have been observed in Antarctic clouds (Lachlan-Cope et al., 2016; O'Shea et al., 2017). Secondary ice processes are believed to magnify ICNCs in polar clouds with important implications for the surface radiative balance (Young et al., 2019), yet the underlying mechanisms remain highly uncertain (Field et al., 2017).

The only well-established secondary ice production (SIP) mechanism that has been extensively implemented in weather forecast and climate models is rime splintering (Hallett and Mossop, 1974), also known as the Hallett–Mossop process (H-M), which refers to the production of ice splinters after collisions of supercooled droplets with ice particles (Hallett and Mossop, 1974; Heymsfield and Mossop, 1984). This process is effective only in a limited temperature range, between −8 and −3 ^∘C, and requires the presence of supercooled liquid droplets both smaller and larger than 13 and 24 µm, respectively (Mossop and Hallett, 1974; Choularton et al., 1980). However, recent studies have shown that H-M cannot sufficiently explain the enhanced ICNCs observed in both Arctic (Sotiropoulou et al., 2020) and Antarctic (Young et al., 2019) clouds. While some Antarctic studies (Vergara-Temprado et al., 2018; Young et al., 2019) suggest that the underestimation of ice multiplication in models might be related to uncertainties in the description of the H-M process, we argue that this is likely driven by the fact that almost no models include other SIP mechanisms.

Another SIP mechanism, identified in recent laboratory studies (Leisner et al., 2014; Lauber et al., 2018), is the generation of ice fragments from the shattering of relatively large frozen drops. This process, however, while very efficient in convective clouds (Korolev et al., 2020), has been found to be ineffective in polar regions (Fu et al., 2019; Sotiropoulou et al., 2020). This is in agreement with Lawson et al. (2017) and Sullivan et al. (2018a) who have shown that drop shattering occurs in clouds with a relatively warm cloud base.

Mechanical breakup (BR) of ice particles that collide with each other is another process that results in ice multiplication (Vardiman, 1978; Takahashi et al., 1995), and it operates over a wide temperature range with maximum efficiency around −15 ^∘C. The limited knowledge of the BR mechanism comes from a few laboratory experiments (Vardiman, 1978; Takahashi et al., 1995) and small-scale modeling (Fridlind et al., 2007; Yano and Phillips, 2011, 2016; Phillips et al, 2017a, b; Sullivan et al. 2018a; Sotiropoulou et al., 2020). To the authors' knowledge only two attempts have been made to incorporate this process into mesoscale models (Sullivan et al., 2018b; Hoarau et al., 2018). Specifically, Hoarau et al. (2018) assumed a constant number of fragments (F_BR) generated per snow–graupel collision in the Meso-NH model, while Sullivan et al. (2018b) implemented a temperature-dependent relationship for F_BR in COSMO-ART for several types of collisions (e.g., crystal–graupel, graupel–hail, etc) based on the results of Takahashi et al. (1995). Phillips et al. (2017a) recently developed a physically based description of F_BR, which is a function of collisional kinetic energy and accounts for the effect of the colliding particles' size and rimed fraction (Ψ). While being more advanced than any other parameterization proposed for BR, this scheme has never been implemented in mesoscale models before; it has only been tested in small-scale models for convective clouds (Phillips et al., 2017b; Qu et al., 2020) and Arctic stratocumulus clouds (Sotiropoulou et al., 2020).

Sotiropoulou et al. (2020) recently showed that the observed ICNCs in Arctic clouds within the H-M temperature zone can be explained only by the combination of BR with the H-M process, which results in a 10- to 20-fold enhancement of the primary ice crystals. Based on their results, we postulate that BR may also play a critical role in Antarctic clouds and can potentially explain the discrepancy between the observed and modeled ICNCs in the region (Young et al., 2019). To test this hypothesis, we implement parameterizations of the BR process in the Morrison microphysics scheme (Morrison et al., 2005) (hereafter M05) in the Weather and Research Forecasting (WRF) model V4.0.1 and examine their influence on the Antarctic clouds observed during the Microphysics of Antarctic Clouds (MAC) field campaign (O'Shea et al., 2017; Young et al., 2019).

2 Observations

2.1 MAC instrumentation

The MAC field campaign was conducted in November–December 2015 over coastal Antarctica and the Weddell Sea with the aim to offer detailed measurements of the microphysical and aerosol properties of the coastal Antarctic atmosphere. MAC included an extensive suite of airborne and ground-based instruments, a detailed description of which can be found in O'Shea et al. (2017). Here we only offer a brief recap of the instrumentation used in this study.

Cloud particle size distributions were derived using the images from a 2D Stereo (2DS; SPEC Inc., USA; Lawson et al., 2006) probe with a nominal size range from 10 to 1280 µm (10 µm pixel resolution). Shattering effects at the probes' inlet were corrected by applying “anti-shatter” tips (Korolev et al., 2011) and inter-arrival time (IAT) post analysis (Crosier et al., 2011). The 2DS is a single particle instrument, measuring all particles that pass through its sample volume, which depends on particle size and the data integration period. For example, at 300 µm, one count measured over a 1 s averaging window equals a concentration of 0.27 L⁻¹; the uncertainty due to counting statistics is 100 %. Total uncertainty is even higher but cannot be quantified.

Aerosol particle measurements of sizes from 0.25 to 32 µm were made using the GRIMM optical particle counter (GRIMM model 1.109), while a cloud aerosol spectrometer (CAS; Baumgardner et al., 2001; Glen and Brooks, 2013) measured particles between 0.6 and 50 µm. Following the methodology of Young et al. (2019) and O'Shea et al. (2017), we only consider GRIMM measurements of particles with diameter smaller than 1.6 µm in our analysis to avoid including data subject to inlet losses at larger particle sizes. Finally, the aircraft also included instrumentation to measure temperature, turbulence, humidity, radiation and surface temperature (King et al., 2008).

2.2 Case study

For our investigations, we focus on the MAC case examined in Young et al. (2019) for which they showed that the H-M process, as currently parameterized in WRF, cannot explain the observed ICNCs. Young et al. (2019) utilized measurements from two MAC flights, M218 and M219, combined in one case study; both flights were conducted on 27 November 2015 over the Weddell Sea (Fig. 1): M218 between 15:18 and 16:42 UTC and M219 between 20:27 and 22:30 UTC. On that day, a low pressure system persisted over the eastern Weddell Sea resulting in a southeasterly flow reaching the aircraft with air mass back trajectories from around the low pressure system, towards the Antarctic Peninsula and southern Patagonia (O'Shea et al., 2017).

Figure 1 Map of Antarctic domains. Colors indicate terrain heights (green to yellow) and sea ice concentrations (blue to white), whereas the purple contours correspond to 500 hPa geopotential heights from the control (CNTRL) simulation at 18:00 UTC, 27 November 2015. The solid black line delimits the 1 km horizontal grid spacing domain, while the dashed one outlines the subset of the nest used for direct comparison with the aircraft data. Orange and blue lines indicate the flight tracks, while the red circle represents Halley station. The small figure in the top left corner indicates the location of the 1 km horizontal grid spacing domain relative to the Antarctic continent.

The temperature and microphysical conditions encountered during these flights are representative of the MAC campaign (see Table 1 in O'Shea et al., 2017, and Fig. S6 in Young et al., 2019). Cloud measurements were collected mainly within the lowest 1.1 km above sea level (a.s.l.) during both flights and within a temperature range of $\sim -\mathrm{9}$ to −3 ^∘C. The sampled stratocumulus clouds were dominated by supercooled liquid drops, while ice formation occurred in isolated ice patches characterized by substantially enhanced ICNCs: the mean (max) ICNCs in these cloud regions were 1.16 (9.03) L⁻¹ and 3.33 (87.31) L⁻¹ for M218 and M219, respectively. The mean concentration of aerosols with sizes between 0.5 and 1.6 µm was 0.56 and 0.41 scm⁻³ (cm⁻³ at standard temperature and pressure) for the two flights. Such low aerosol conditions and concurrently high ICNC concentrations within this temperature range are frequently found in the western Antarctic Peninsula (Lachlan-Cope et al., 2016). Moreover, similar cloud droplet concentrations (N_drop) were measured during both flights (Young et al., 2019): the mean N_drop was 82.7 cm⁻³ for M218 and 100.4 cm⁻³ for M219, which are comparable with previous observations from the Antarctic Peninsula (Lachlan-Cope et al., 2016).

3 Modeling methods

3.1 Model setup

This study is conducted with the WRF model (Skamarock et al., 2008) version 4.0.1 by applying the same model setup as in Young et al. (2019). Two domains with a respective horizontal resolution of 5 and 1 km are used, in which the inner one is two-way nested to the parent domain (Fig. 1). The polar stereographic projection is applied. The outer domain is centered at 74.2^∘ N, 30^∘ E and includes 201 × 201 grid points, while the second domain consists of 326 × 406 grids. Both domains have a high vertical resolution with 70 eta levels, 25 of which correspond to the lowest 2 km of the atmosphere. The model top is set to 50 hPa. The simulation period spans from 26 to 28 November 2014, 00:00 UTC, allowing for a 24 h spin up period before the day of interest (27 November). The model time step is set to 30 (6) s for the outer (inner) domain, while output data are produced every 30 min.

Input data for the initial, lateral and boundary conditions for the simulations are obtained from the European Centre for Medium-Range Weather Forecasting reanalysis (Dee et al., 2011), as recommended by Bromwich et al. (2013). For both shortwave and longwave radiation components, the RRTMG radiation scheme (Rapid Radiative Transfer Model for general circulation models) is applied. The Mellor–Yamada–Nakanishi–Niino (MYNN; Nakanishi and Niino, 2006) level-2.5 closure planetary boundary layer (PBL) and surface options are also implemented in combination with the Noah Land Surface Model (Noah LSM; Chen and Dudhia, 2001) which includes a simplified thermodynamic sea ice model. Given the short run length, time-varying sea ice concentrations are not utilized. Young et al. (2019) used the Polar WRF V3.6.1 to represent fractional sea ice, a capability not available in standard WRF V3.6. However, this option has been made available in the more recent V4.0.1 that we use in this study. Following Young et al. (2019), the sea ice albedo is set to 0.82 with a default thickness of 3 m, and snow accumulation depth on sea ice is allowed to vary between 0.001 and 1.0 m.

A so-called “cumulus parameterization” for shallow-convection subgrid processes is not activated in both domains to ensure all cloud processes are represented by the grid-scale microphysics scheme. Note that 5 km is a general upper limit for a convection-resolving resolution (Klemp, 2006; Prein et al., 2015). Cloud microphysics are parameterized following Morrison et al. (2005), hereafter M05. M05 performs well in reproducing Antarctic clouds, resulting in improved representation of the liquid phase and thus the cloud radiative effects being compared to less advanced microphysical schemes (Listowski and Lachlan-Cope, 2017; Hines et al., 2019). This bulk microphysics scheme predicts mixing ratios and number concentrations for cloud ice, rain, snow and graupel species. While the mass mixing ratio of cloud water is a prognostic variable, N_drop is a constant parameter. The default value of the scheme is 200 cm⁻³; here N_drop is set to 92 cm⁻³, which is the mean value of M218 and M219 flight measurements (see Sect. 2.2).

3.2 Sensitivity simulations

A detailed description of the ice formation processes in M05 and the implemented BR parameterizations is offered in Appendix A and B, respectively. We assume that collisions that include at least one large particle (thus ice–snow, ice–graupel and graupel–snow, snow–snow, and graupel–graupel) result in ice multiplication; contributions from collisions between small ice particles (cloud ice) are neglected. In addition to the control (CNTRL) simulation, which corresponds to the default setup of M05 and accounts only for H-M, we perform seven sensitivity simulations with varying descriptions of F_BR. We also perform an additional simulation as in CNTRL except with no H-M and thus no SIP at all, which is referred to as NOSIP in the text.

In two sensitivity simulations with active breakup, we assume, as in Hoarau et al. (2018), a constant number of fragments generated per collision. This number is constrained by in situ measurements from the Arctic (Schwarzenboeck et al., 2009) which indicated that one-branch ice crystals are more common in polar clouds, resulting in the ejection of a single fragment after collision with another ice particle. However, this analysis (Schwarzenboeck et al., 2009) included only dendritic crystals with a size larger than 300 µm. Based on these results, we perform two simulations: FRAG1 assumes all collision types generate one fragment without any size restrictions, while FRAG1siz allows for ice multiplication only if the particle that undergoes fragmentation is larger than 300 µm. Note that because cloud ice with a characteristic diameter larger than 250 µm is converted to snow in the M05 scheme, collisions that include cloud ice are assumed to not result in any multiplication in FRAG1siz.

The standard temperature-dependent formula of Takahashi et al. (1995) for F_BR, applied in Sullivan et al. (2018b), is tested here in the TAKAH simulation. However, Takahashi et al. (1995) used 2 cm hail balls in their experiments, which is an unrealistic setup. For this reason, we perform an additional simulation, TAKAHsc, in which this relationship is further scaled with size (see Appendix B).

Finally, the Phillips parameterization is implemented in three simulations with varying rimed fractions (Ψ) for the cloud ice or snow particles that undergo fragmentation; Ψ is not predicted in most bulk microphysics schemes, including M05, and thus it is prescribed as a constant. Note that F_BR is a function of Ψ only for the ice crystals or snowflakes that undergo breakup but not for graupel (Appendix B). Graupel is assumed to have Ψ≥0.5, while the other ice types are characterized by a lower rimed fraction. For this reason, we will consider values of Ψ for cloud ice and snow between 0.2 (lightly rimed) and 0.4 (heavily rimed) (Phillips et al., 2017a, b). These simulations are referred to as PHIL0.2, PHIL0.3 and PHIL0.4 in the text, in which the number indicates the assumed values of Ψ for cloud ice and snow.

Figure 2 (a) Mean ice number concentrations (cloud ice + snow + graupel; N_isg) as a function of temperature for the whole MAC campaign (pink), our case study (gray) and the eight model simulations. Gray dots indicate point observations. (b) Relative frequency distribution of N_isg, binned in 0.5 L⁻¹ intervals and scaled with maximum frequency. Ice properties are calculated for particles > 80 µm and for N_isg > 0.005 L⁻¹ within the lowest 1.5 km a.s.l.

4 Results

4.1 BR effect on microphysical properties

In Fig. 2a, the modeled total ice number concentrations (cloud ice + snow + graupel; N_isg) derived for the region encompassing the two MAC flights (Fig. 1) are compared with measurements derived from the 2D Stereo (2DS) probe (see Sect. 2.1 for details). ICNCs in Fig. 2 are interpolated to match the time resolution of the temperature measurements. Then cloud ice statistics are calculated for N_isg > 0.005 L⁻¹, an indicator for the presence of an ice patch (O'Shea et al., 2017; Young et al., 2019). Moreover, since 2DS cannot resolve the shape (thus cloud phase) of particles smaller than 80 µm, only modeled ice particles with sizes larger than this threshold are considered in Fig. 2, like in Young et al. (2019). While mean and maximum statistics are discussed below, additional statistical metrics (e.g., median and interquartile range) are shown in Fig. S1 (Sect. S1).

The mean observed N_isg for the whole MAC campaign generally fluctuates between 0.5 and 4.5 L⁻¹. The variation in N_isg with temperature is somewhat larger for our case study (27 November) as maximum mean concentration goes up to ∼ 6.4 L⁻¹ at $T=-$6.5 ^∘C. Consistently lower concentrations are observed for temperatures ≤ −7 ^∘C, but the temperature statistics are poor for this temperature range as very few observations are available (Fig. 2a). The CNTRL simulation consistently underestimates the mean observations, producing mean N_isg∼0.1 L⁻¹ over the examined temperature range (Fig. 2a). NOSIP produces similar results to CNTRL, suggesting that the H-M process included in default M05 (CNTRL) is hardly effective at all.

PHIL0.2 and PHIL0.3 also produce similar mean ICNCs to CNTRL (Fig. 2a, b), suggesting that lightly to moderately rimed ice particles do not contribute to ice multiplication through collisional breakup. However, the higher rimed fraction in PHIL0.4 results in very good agreement with mean observations (Fig. 2a). FRAG1siz produces very weak ice multiplication that cannot explain observed ICNCs, but when the size restrictions are ignored (FRAG1), the model usually overestimates mean values, particularly at temperatures ≤ 7 ^∘C for which no good measurement statistics are available (see discussion above). A similar behavior is shown by TAKAHsc, while TAKAH is the simulation that produces the most unrealistically high mean ICNCs (Fig. 2a).

Overall, CNTRL, PHIL0.2 and PHIL0.3 cannot reproduce the observed spectrum (Fig. 2b) and substantially underestimate the frequency of ICNCs larger than 2 L⁻¹. FRAG1siz gives a wider representation of the size range in better agreement with observations; however, it cannot reproduce ICNCs larger than 30 L⁻¹. PHIL0.4, FRAG1 and TAKAHsc can successfully reproduce the observed range of values (Fig. 2b), but the relative frequency remains substantially more underestimated in PHIL0.4. In contrast, TAKAH often overestimates the relative frequency of ICNC values larger than 10 L⁻¹ (Fig. 2b). Maximum ICNCs in TAKAH are 2952 L⁻¹, which is about 34 times larger that the observed maximum value: 88 L⁻¹. The maximum ICNC in PHIL0.4 is 161 L⁻¹, which agrees to within a factor of 2 with observations, while it is slightly larger for FRAG1 and TAKAHsc: 188 and 181 L⁻¹, respectively. On the other hand, maximum concentrations are substantially underestimated in CNTRL (5.8 L⁻¹), PHIL0.2 (2 L⁻¹) and PHIL0.3 (3.5 L⁻¹).

Vertical distributions of cloud ice (N_i), graupel (N_g), snow (N_s) and total ICNC (N_isg) number concentrations are examined in Fig. 3a–d for all simulations except TAKAH, which produces unrealistically large concentrations. The observed ICNCs are also shown in Fig. 3a and c. For consistency with M05, which converts all cloud ice particles with characteristic diameters larger than 250 µm to snow, the same threshold is adopted for splitting the observational dataset in these two ice categories. Graupel concentrations cannot be distinguished in the measurements (hence no “Nov 27” profile in Fig. 3b); however, the model simulations that are in better agreement with observations (Fig. 2) suggest that these are negligible compared to cloud ice and/or snow concentrations. Graupel concentrations in Fig. 3b are shown for the whole size spectrum. In contrast, cloud ice (Fig. 3a), snow (Fig. 3c) and total ICNCs (Fig. 3d) include only particles with a size larger than 80 µm for consistency with the observations shown in the same panel.

Figure 3 Mean vertical profiles of number concentrations of modeled (a) cloud ice, (b) graupel, (c) snow and (d) total ICNCs for seven simulations. Gray lines represent measured concentrations with diameters (a) smaller and (c) larger than 250 µm. Graupel concentrations cannot be distinguished in the measurements (hence no gray profile in panel b). Ice properties from the model are calculated for N_isg > 0.005 L⁻¹. For consistency with observations, only particles with sizes > 80 µm are included in the modeled profiles in panels (a), (c) and (d).

PHIL0.2 and PHIL0.3 produce similar N_i (Fig. 3a) to CNTRL but reduced N_g (Fig. 3b) values and reduced N_s (Fig. 3c); these mean N_i and N_s profiles are orders of magnitude lower than the observed values. N_i (Fig. 3a) and total ICNCs (Fig. 3d) are enhanced in FRAG1siz compared to CNTRL but remain a substantially lower than the observed values. PHIL0.4 produces N_i close to the observations (Fig. 3a), while N_s (Fig. 3c) and thus total ICNCs (Fig. 3d) are somewhat underestimated. FRAG1 and TAKAHsc are in better agreement with N_s observations (Fig. 3c), especially at heights below 750 m a.s.l; this is also reflected in total ICNCs (Fig. 3d). Activating BR generally results in a reduction of N_g (Fig. 3b). This decrease can be larger than 1 order of magnitude in some of the best performing simulations compared to CNTRL; however, we cannot assess which of these graupel profiles better represents reality. Nevertheless, we can overall conclude that PHIL0.4, FRAG1 and TAKAHsc result in improved agreement of the vertical distribution of total ICNCs with observations compared to the rest of the simulations (Fig. 3d), including the default setup of M05. Moreover, cloud ice concentrations (Fig. 3a) are comparable to snow concentrations (Fig. 3c) in these three simulations, in agreement with observations. In contrast, simulations with deactivated or negligible BR result in substantially larger number of snow than cloud ice particles. This indicates that BR shifts the ice particle spectra to smaller sizes, which results in a more realistic representation of the ice microphysical characteristics.

Figure 4 (a) Mean liquid water content (LWC) as a function of temperature for our case study (gray) and the eight model simulations. (b) Relative frequency distribution of LWC, binned in 0.01 g m⁻³ intervals, scaled with maximum frequency. Only values greater than 0.01 g m⁻³ within the lowest 1.5 km a.s.l. are included in the analysis.

The simulated liquid water content (LWC) is compared with CAS observations in Fig. 4. All simulations, except TAKAH, produce similar or slightly overestimated mean LWC at temperatures $\le -\mathrm{3.5}$ ^∘C; at −3 ^∘C, the mean observed values are higher (Fig. 4a). An overestimation of LWC in these runs is more evident in Fig. 4b; the observed spectrum does not include values larger than 0.5 g m⁻², while the simulated spectra are wider. An exception to this is the TAKAH simulation which underestimates mean LWCs and glaciates clouds at temperatures below −7 ^∘C (Fig. 4a), while it produces a narrower LWC spectrum compared to the one observed (Fig. 4b). Apart from TAKAH, the remaining simulations produce similar liquid water properties with minor improvements in the runs with reduced LWC values, e.g., in FRAG1 (Fig. 4b).

4.2 BR effect on surface radiation

To examine how deviations in ICNCs affect climate, mean radiative fluxes at the surface and at the top of the atmosphere (TOA) for all model simulations are presented in Table 1. Note that mass mixing ratio fields for all cloud species are provided from the microphysics to the radiation scheme, but no information on droplet and ice effective radius is exchanged.

Table 1 Mean modeled downward and upward shortwave (SWD_SFC, SWU_SFC) and longwave (LWD_SFC, LWU_SFC) surface radiation, along with upward shortwave and longwave (SWU_TOA, LWU_TOA) radiation at the top of the atmosphere, during flights M218 and M219. Model results are averaged over the dashed rectangular area in Fig. 1.

Increasing BR multiplication has a pronounced impact on shortwave radiation as it results in decreasing sunlight reflection and thus increasing downward surface radiation (SWD_SFC). Upward surface radiation (SWU_SFC) is a function of SWD and thus exhibits similar behavior. This is due to the fact that increased BR efficiency (Fig. 2) results in decreased liquid water path (LWP) and cloud albedo. The difference between CNTRL and the simulations that improve ICNC representation (PHIL0.4, FRAG1 and TAKAHsc) fluctuates between 12.5 and 24.4 W m⁻² for SWD_SFC and 6.9 and 13.7 W m⁻² for SWU_SFC (Table 1).

Cloud longwave radiative effects are mainly determined by cloud liquid properties since liquid water is more opaque to longwave radiation than ice particles. However, no substantial differences in mean LWP are indicated for CNTRL (40.1 g m⁻²), PHIL0.2 (45.6 g m⁻²), PHIL0.3 (46.2 g m⁻²), PHIL0.4 (32.5 g m⁻²) and FRAG1siz (38.6 g m⁻²) since LWP values fall within the black body emission range (Stephens, 1978). Optically thinner clouds are produced in TAKAHsc (20.7 g m⁻²) and especially in FRAG1 (21.1 g m⁻²) and TAKAH (3.2 g m⁻²) runs. Note that most simulations, including CNTRL, produce wider LWC spectra than those observed by overestimating cloud liquid (Fig. 4b). Generally, decreasing liquid content is in better agreement with observations (see Sect. 4.1), suggesting that including the BR process in M05 likely shifts the simulated LWPs towards more realistic values. However, excessive ice multiplication, as in TAKAH, results in unrealistic liquid properties (Fig. 4a) and thus errors in surface radiation.

Pronounced reduction in LWD_SFC is only found for the simulations FRAG1, TAKAH and TAKAHsc, which have a mean LWP well below 30 g m⁻², the lowest limit of the black body emission range (Stephens, 1978). In all other simulations, the reduction in cloud liquid due to BR is not large enough to alter the cloud emissivity significantly. The upward longwave component (LWU_SFC) remains almost unaffected in all simulations (≤0.6 W m⁻²). Young et al. (2019) showed that underestimation of ICNCs results in significant positive and negative biases in the surface cloud radiative forcing (CRF) over the coastal areas; our results agree with their findings as CRF biases vary between −68 and +87 W m⁻² for the most realistic simulations (Fig. S2, Sect. S2). Furthermore, the difference between CNTRL and the realistic simulations in upward radiation flux at TOA (Table 1) is also more pronounced for the shortwave component (SWU_TOA), fluctuating between 4 and 8.6 W m⁻², and less significant for LWU_TOA (1.3 and 2.7 W m⁻²). Ultimately, both surface and TOA radiation results indicate that a correct representation of SIP in the atmospheric models is critical for the projection of the energy budget and thus for the future Antarctic climate.

4.3 Sensitivity to uncertainties in H-M description

To investigate the interactions between BR and H-M, we compare simulations in which the H-M efficiency is either enhanced or turned off. Young et al. (2019) remove all liquid thresholds from the H-M description, allowing for the process to become active over the whole droplet spectrum. However, this change resulted in very weak ICNC enhancement in their simulations. Here we further remove all graupel and/or snow thresholds from the H-M description (Appendix A) following the methodology of Atlas et al. (2020), which implies that there are no size restrictions for the initiation of the process. This modification is applied to the CNTRL, PHIL0.3 and PHIL0.4 setups, resulting in three additional sensitivity tests: CNTRL_NOTHRES, PHIL0.3_NOTHRES and PHIL0.4_NOTHRES, respectively. Furthermore, in addition to NOSIP which corresponds to the CNTRL setup but without H-M, another two simulations are performed with BR active but H-M completely deactivated: PHIL0.3_NOHM and PHIL0.4_NOHM.

Figure 5 Total ice number concentrations (N_isg) for particles > 80 µm as a function of temperature for the whole MAC campaign (pink), our case study (gray) and the nine sensitivity simulations with varying treatments of the H-M process (see Sect. 4.3). Mean values and the 95th percentile are shown in panels (a) and (b), respectively.

Mean ICNCs in CNTRL_NOTHRES are enhanced by on average a factor of 3 compared to CNTRL (Fig. 5a). However, this simulation underestimates concentrations at temperatures larger than −7 ^∘C; the mean observed value at this range is 2 L⁻¹, while the simulated mean is 0.3 L⁻¹. Good agreement between CNTRL_NOTHRES and observations is only achieved at temperatures <−7 ^∘C, in which statistical metrics for the two MAC cases are poor (see Sect. 4.1). While PHIL0.3 did not result in any substantial multiplication, mean ICNCs in PHIL0.3_NOTHRES are about 5 times larger. A similar enhancement is observed in PHIL0.4_NOTHRES compared to PHIL0.4, but it becomes substantially larger at colder temperatures. However, while the 95th percentiles for CNTRL_NOTHRES and PHIL0.3_NOTHRES are similar and more comparable to observations, PHIL0.4_NOTHRES produces values larger than 10 L⁻¹ at all temperatures considered (Fig. 5b).

Excluding the temperature range (< −7 ^∘C) that does not include a substantial number of measurements to evaluate model results, mean ICNC observations generally lay between PHIL0.3_NOTHRES and PHIL0.4_NOTHRES in this set of simulations, while CNTRL_NOTHRES produces somewhat lower values (Fig. 5a). However, this setup overestimates H-M efficiency as it does not include any size limitations, which is not consistent with current knowledge on the H-M mechanism derived from laboratory studies (Hallet and Mossop, 1974; Choularton et al., 1980). Nevertheless, the adapted thresholds are ad hoc and unsuitable for the examined conditions (Atlas et al., 2020); these should be refined to get a more a realistic H-M effect in polar clouds.

Deactivating H-M completely does not substantially impact the results. This further confirms the fact that the prescribed ad hoc thresholds prevent the initiation of the process in the studied conditions. Furthermore, it indicates that the BR mechanism can explain the observed ICNCs independently of whether H-M is active or not.

4.4 Sensitivity to uncertainties in primary ice formation

None of the utilized primary ice nucleation parameterizations are likely representative of the pristine conditions encountered over the high-latitude Southern Ocean; thus, it is likely that primary ice formation is overestimated in this case. To examine how the uncertainty in parameterizations for primary ice affects SIP efficiency, we perform two sets of simulations by dividing or multiplying the efficiency of all primary ice production mechanisms (immersion freezing, contact freezing, and deposition and/or condensation-freezing nucleation) by a factor of 10. Specifically, the first set with diminished ice nucleation includes CNTRL_INP0.1, PHIL0.3_INP0.1 and PHIL0.4_INP0.1, while the second set with enhanced nucleation consists of CNTRL_INP10, PHIL0.3_INP10 and PHIL0.4_INP10.

Decreasing primary ice production by a factor of 10 results in a substantial decrease in ICNCs: about 25 %–40 % fewer values (depending on the BR setup) exceed the 0.005 L⁻¹ threshold and are included in the presented statistics (Fig. 6) compared to the simulations with the standard primary ice formation setup. BR multiplication is again evident only in the simulation with a high Ψ (Fig. 6). The ICNC enhancement in PHIL0.4_INP0.1 compared to CNTRL_INP0.1 is about on average 7 times larger (Fig. 6a). While these results are in good agreement with observations in the temperature range with poor measurement statistics (see discussion above), PHIL0.4_INP0.1 underestimates mean observations at temperatures > −7 ^∘C (Fig. 6a). Moreover, the 95th percentile in this simulation remains substantially lower compared to observed values at all temperatures (Fig. 6b). These results indicate that an initial ICNC concentration of 0.1 L⁻¹ is essential to initiate significant multiplication that can reproduce the observed values.

Figure 6 Same as in Fig. 5 but for sensitivity simulations with varying INP conditions (see Sect. 4.4).

Increasing primary ice production by an order of magnitude in CNTRL_INP10 still results in underestimated ice concentrations than those observed, providing additional evidence for the significant role of SIP in these conditions. The increased concentration of primary ice crystals enhances BR efficiency in PHIL0.3_INP10 compared to PHIL0.3; however, the produced mean concentrations still are lower than the observed. The difference in mean ICNCs between PHIL0.4_INP10 and PHIL0.4 is on average about a factor of 2; thus, PHIL0.4_INP10 is also close to the observations, especially those from the whole MAC campaign (Fig. 6a). N_isg95 is also larger in PHIL0.4_INP10 compared to PHIL0.4 and in better agreement with observations (Fig. 6b).

In summary, the above results indicate that BR is not sufficient to explain observations when the available primary ice concentrations are substantially lower than 0.1 L⁻¹, which is the mean primary ICNC produced in the NOSIP simulation (Fig. 2a). Yet, INPs over the Southern Ocean are often substantially lower (McCluskey et al., 2018; Schmale et al., 2019; Welti et al., 2020). Ice seeding from clouds above the boundary layer was suggested by Young et al. (2019) as a key contributor to the primary ICNC levels for the studied case (see their Supplement). Another process that can enhance primary ice nucleation is aerosol transport from the Antarctic continent, where terrestrial INPs are expected to be higher (Vergara-Temprado et al., 2018). Moreover, a combination of BR and the H-M mechanism, whose efficiency is substantially restricted in the current version of M05 (see Sect. 4.3), might also explain the observed number concentrations; this was also the case for Arctic stratocumulus conditions in Sotiropoulou et al. (2020). Understanding these interactions between different processes in the Antarctic region would likely provide insights into the conditions that favor the development of isolated ice patches with substantially high ICNCs within predominantly supercooled liquid clouds. In higher INP conditions, which are likely less representative of the coastal Antarctic climate, the activation of BR parameterization still results in improved model performance.

5 Conclusions

Our results indicate that collisional breakup of ice crystals can explain observations of enhanced ICNCs in coastal Antarctic clouds as long as ∼ 0.1 L⁻¹ of primary ice crystals are available (as produced in NOSIP simulation). This likely is a key threshold that can lead the development of isolated ice patches with substantially high ICNCs in predominantly supercooled liquid clouds (Grosvenor et al., 2012; O'Shea et al., 2017). Over the Southern Ocean, when INPs are generally sparse (McCluskey et al., 2018; Schmale et al., 2019; Welti et al., 2020), such conditions could likely be achieved through ice seeding (as likely happens in the examined case) or through INP transport from the Antarctic continent where INP concentrations are predicted to be higher (Vergara-Temprado et al., 2018). These points remain for future confirmation.

Although BR has been observed in polar conditions before (Rangno and Hobbs, 2001; Schwarzenboeck et al., 2009), this mechanism is currently not implemented in most weather prediction and climate models. The more advanced Phillips et al. (2017a) parameterization produces realistic ICNCs in Antarctic clouds as long as a high rimed fraction is prescribed for the particles that undergo fracture, which is in agreement with Sotiropoulou et al. (2020). This indicates that our conclusions may not hold for winter clouds in the region which contain less supercooled liquid water (Listowski et al., 2019) and are less prone to riming. However, for the studied case, a comparison of vapor deposition rates with riming rates (which include mass changes due to collisions with droplets and/or raindrops and due to contact and/or immersion freezing) for the CNTRL simulation indicates that these two are on average comparable for cloud ice, while riming rates are substantially larger than vapor deposition rates for snow (not shown). These results suggest that prescribing a high rimed fraction for cloud ice and snow in M05 is not unreasonable; nevertheless Ψ in reality is highly variable for different temperature and microphysical conditions. The more simplified parameterization by Takahashi et al. (1995) also produces improved results as long as the dependence of F_BR on the ice particle size is accounted for. The results of this setup are similar to the parameterization that assumes that all ice particle collisions generate one fragment.

The very few existing BR descriptions in mesoscale models either do not account for size limitations (Sullivan et al., 2018b) or do not account for all collision types (Hoarau et al., 2018), which limits their realism. Increasing ICNCs from BR alters significantly the radiative effects of summer mixed-phase Antarctic clouds; these clouds play a critical role in the surface melting of ice shelves in the vicinity of the Weddell Sea (Gilbert et al., 2020), and thus their accurate microphysical representation in models is of great importance.

Appendix A: Ice formation processes in M05 scheme

The standard M05 scheme includes three primary ice production mechanisms through heterogeneous nucleation (immersion freezing, contact freezing, and deposition and/or condensation-freezing nucleation) and one SIP process (H-M).

Immersion freezing of cloud droplets and rain is based on the work of Bigg (1953). This mechanism is active below −4 ^∘C and produces a raindrop freezing rate that depends on the degree of supercooling and the number concentration and volume of supercooled drops. The Meyers et al. (1992) description is used for contact freezing also active below −4 ^∘C. The effective diffusivity of the contact nuclei to the drops are estimated from Brownian motion similar to Young (1974): $D_{\mathrm{ap}}=RT\left(\mathrm{6}p{r}_{\mathrm{i}}{N}_{\mathrm{A}}m{)}^{-\mathrm{1}}\right[\mathrm{1}+\mathrm{0.0737}T(\mathrm{2880}P{)}^{-\mathrm{1}}{r}_{\mathrm{i}}^{-\mathrm{1}}$], where R is the universal gas constant, N_A is Avogadro's number, m is the dynamic viscosity of air, T is the temperature and P is the air pressure, and the radius of ice nuclei r_i is assumed to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ m. The factor in the square brackets is a correction factor accounting for the mean free path of air molecules relative to the size of the ice nuclei (all units are in the meter, kilogram and second, MKS, system of units).

The default parameterization for deposition and/or condensation-freezing ice nucleation in M05 is from Cooper (1986), which depends only on temperature and is active below −8 ^∘C in liquid saturated conditions or when ice supersaturation exceeds 8 %. However, the aerosol-aware DeMott et al. (2010) parameterization for heterogeneous nucleation has been shown to compare better with in-cloud ice measurements over the Antarctic Peninsula than Cooper (Listowski and Lachlan-Cope, 2017), although none of these schemes is likely accurate as they have not been developed for such pristine conditions. Nevertheless, both Cooper and DeMott parameterizations produce similar primary ice concentrations over the temperature range covered by the observations, but Cooper predicts more primary ice at lower temperatures (< 13 ^∘C), which might affect the representation of higher-altitude clouds (see Supplement in Young et al., 2019). For this reason, we apply the DeMott description in our simulations, in which the mean aerosol concentration of particles with sizes between 0.5 and 1.6 µm for the two flights (0.49 scm⁻³) is used as input (Young et al., 2019). Uncertainty in this formulation is investigated through a number of sensitivity tests (Sect. 4.4).

The H-M parameterization, adapted from Reisner et al. (1998), is based on the laboratory experiments conducted by Hallett and Mossop (1974), who found a maximum of ∼ 350 splinters per milligram of rime generated at around −5 ^∘C:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{{\mathrm{dNi}}_{\mathrm{HM}}}{\mathrm{d}t}}=\mathit{\rho }{\mathrm{3.510}}^{\mathrm{8}}f\left(T\right){\displaystyle \frac{{\mathrm{dm}}_{\mathrm{rime}}}{\mathrm{d}t}},\end{array}$$

where dNi_HM ∕ dt is the number of new fragments produced at a given time step, f(T) is the temperature-dependent efficiency of the process, ρ is the air density and dm_rime ∕ dt is the mass production rate of rime on snow or graupel due to the accretion of cloud and rain drops; f(T) is 0 for T < −8 ^∘C and T > −3 ^∘C, 1 for T= −5 ^∘C, increasing linearly between these two extremes for $T\ge -\mathrm{8}$ ^∘C and $T\le -\mathrm{3}$ ^∘C.

Furthermore, for H-M to become activated in M05, two conditions must be met: (a) snow (or graupel) mass mixing ratios must be greater than 0.1 g kg⁻¹ and (b) cloud liquid (or rain) water mass mixing ratios should be greater than 0.5 (or 0.1) g kg${}^{-\mathrm{1}}.$To achieve a good agreement between modeled and observed ICNCs for the simulated case, Young et al. (2019) had to remove condition (b) and multiply the H-M efficiency by a factor of 10.

Appendix B: Parameterizing collisional breakup in M05

There are three types of ice particles considered in the M05 scheme: small (cloud) ice, snow and graupel. Ice multiplication is allowed after cloud ice–snow, cloud ice–graupel, graupel–snow, snow–snow and graupel–graupel collisions. The standard M05 scheme includes a description for collisions between cloud ice and snow to represent the accretion process following the “continuous collection” approach:

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{AC}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\mathit{\rho }{E}_{\mathrm{col}}\mathrm{\Gamma }\left(b_{\mathrm{s}}+\mathrm{3}\right)a_{\mathrm{s}}{\displaystyle \frac{N_i{\mathrm{N}}_{\mathrm{0}s}}{{\mathit{\lambda }}_{\mathrm{s}}^{\left(b_{\mathrm{s}}+\mathrm{3}\right)}}},\end{array}$$

where dN_iAC ∕ dt is the rate of ice crystal number concentration collected by snow, N_0s and λ_s are the intercept and slope parameters of the snow size distribution represented by an inverse exponential function, Γ is the Euler gamma function, and a_s and b_s are the characteristic parameters for snow in the fall speed–diameter relationship (Morrison et al., 2005); a_s includes a density correction factor (Heymsfield et al., 2007). Note that the diameter (d_i) and terminal velocity (u_i) of cloud ice particles are considered much smaller than those of snow: d_i≪d_s and u_i≪u_s, so that they are neglected in Eq. (2). E_col is the collection (sticking) efficiency between ice particles, set to 0.1; hence, it is assumed that only 10 % of cloud ice particles that collide with snow are actually collected. We assume the remaining 90 % of collisions result in ice particle breakup; hence, the following relationship gives the rate of cloud ice–snow collisions that contribute to ice multiplication:

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{is}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\mathit{\rho }(\mathrm{1}-E_{\mathrm{col}})\mathrm{\Gamma }\left(b_{\mathrm{s}}+\mathrm{3}\right)a_{\mathrm{s}}{\displaystyle \frac{N_i{\mathrm{N}}_{\mathrm{0}s}}{{\mathit{\lambda }}_{\mathrm{s}}^{\left(b_{\mathrm{s}}+\mathrm{3}\right)}}}.\end{array}$$

In the default M05, collisions between cloud ice and graupel particles are neglected as it is assumed that the collection efficiency of such collisions is negligible. To represent cloud ice–graupel collisions for ice multiplication, we use Eq. (B2), but the size distribution and fall speed parameters of snow are replaced by those for graupel. Moreover, since cloud ice is not collected by graupel particles, we assume that 100 % of these collisions result in cloud ice breakup:

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{ig}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\mathit{\rho }\mathrm{\Gamma }\left(b_{\mathrm{g}}+\mathrm{3}\right)a_{\mathrm{g}}{\displaystyle \frac{N_i{N}_{\mathrm{0}g}}{{\mathit{\lambda }}_{\mathrm{g}}^{\left(b_{\mathrm{g}}+\mathrm{3}\right)}}}.\end{array}$$

In the default M05 scheme, collisions between snow and graupel are also neglected because it is assumed that the collection efficiency for such collisions is negligible. For this study, graupel–snow collisions are treated using expressions similar to those for raindrop–snow collisions in M05. These are adapted from Ikawa and Saito (1991) and represent collisions between relatively large precipitation-sized particles:

$$\begin{array}{}\text{(B4)}& {\displaystyle \frac{\mathrm{d}{Q}_{i_{\mathrm{sg}}}}{\mathrm{d}t}}={\mathit{\pi }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{s}}\mathit{\rho }\left|\mathrm{\Delta }{u}_{m_{\mathrm{sg}}}\right|{\displaystyle \frac{N_{\mathrm{0}s}{N}_{\mathrm{0}g}}{{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{3}}}}\left({\displaystyle \frac{\mathrm{5}}{{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{3}}{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{2}}{{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{0.5}}{{\mathit{\lambda }}_{\mathrm{s}}{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{3}}}}\right),\text{(B5)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{sg}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\mathit{\rho }\left|\mathrm{\Delta }{u}_{n_{\mathrm{sg}}}\right|N_{\mathrm{0}s}{N}_{\mathrm{0}g}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{3}}{\mathit{\lambda }}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{s}}{\mathit{\lambda }}_{\mathrm{g}}^{\mathrm{3}}}}\right),\end{array}$$

where

$$\begin{array}{}\text{(B6)}& \left|\mathrm{\Delta }{u}_{m_{\mathrm{sg}}}\right|={\left({\left(\mathrm{1.2}{u}_{\mathrm{ms}}-\mathrm{0.95}{u}_{\mathrm{mg}}\right)}^{\mathrm{2}}+\mathrm{0.08}{u}_{\mathrm{mg}}{u}_{ms}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(B7)}& \left|\mathrm{\Delta }{u}_{n_{\mathrm{sg}}}\right|={\left(\mathrm{1.7}{\left(u_{\mathrm{ns}}-u_{\mathrm{ng}}\right)}^{\mathrm{2}}+\mathrm{0.3}{u}_{\mathrm{ng}}{u}_{\mathrm{ns}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where dQ_isg ∕ dt and dN_isg ∕ dt represent the bulk rates at which snow mass and number concentration collide with graupel and contribute to ice multiplication through fragmentation. Corrections in the mass-weighted (or number-weighted) difference in terminal velocity $\mathrm{\Delta }{u}_{m_{\mathrm{sg}}}\left(\mathrm{or}\mathrm{\Delta }{u}_{n_{\mathrm{sg}}}\right)$ of the colliding particles (Eqs. B6, B7) are adapted from Mizuno (1990) and Reisner et al. (1998) to account for underestimates when u_ns≈u_ng.

M05 also includes a description for collisions between snowflakes to represent snow aggregation that is based on Passarelli (1978):

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{\mathrm{d}{N}_{{\mathrm{s}}_{\mathrm{AG}}}}{\mathrm{d}t}}={\displaystyle \frac{-\mathrm{1108}{a}_{\mathrm{s}}{E}_{\mathrm{col}}}{\mathrm{4}\times \mathrm{720}}}{\mathit{\pi }}^{\frac{\mathrm{1}-b_{\mathrm{s}}}{\mathrm{3}}}{\mathit{\rho }}^{\frac{\mathrm{2}+b_{\mathrm{s}}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{s}}^{\frac{-\mathrm{2}-b_{\mathrm{s}}}{\mathrm{3}}}{Q}_{\mathrm{s}}^{\frac{\mathrm{2}+b_{\mathrm{s}}}{\mathrm{3}}}{N}_{\mathrm{s}}^{\frac{\mathrm{4}-b_{\mathrm{s}}}{\mathrm{3}}}.\end{array}$$

Based on this expression we parameterize the number of snow–snow collisions that contribute to ice multiplication as follows:

$$\begin{array}{}\text{(B9)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{ss}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1108}{a}_{\mathrm{s}}(\mathrm{1}-E_{\mathrm{col}})}{\mathrm{4}\times \mathrm{720}}}{\mathit{\pi }}^{\frac{\mathrm{1}-b_{\mathrm{s}}}{\mathrm{3}}}{\mathit{\rho }}^{\frac{\mathrm{2}+b_{\mathrm{s}}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{s}}^{\frac{-\mathrm{2}-b_{\mathrm{s}}}{\mathrm{3}}}{Q}_{\mathrm{s}}^{\frac{\mathrm{2}+b_{\mathrm{s}}}{\mathrm{3}}}{N}_{\mathrm{s}}^{\frac{\mathrm{4}-b_{\mathrm{s}}}{\mathrm{3}}}.\end{array}$$

Because snow aggregation does not result in any mass transfer, the snow mass involved in these collisions is not calculated by the default M05 scheme. We obtain a description of dQ_iss ∕ dt by applying the size distribution and fall speed parameters of snow in the analytical solution for self-collection derived by Verlinde et al. (1990):

$$\begin{array}{}\text{(B10)}& {\displaystyle \frac{\mathrm{d}{Q}_{i_{\mathrm{ss}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{914}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{48}\mathit{\rho }{\mathit{\rho }}_{\mathrm{s}}}}(\mathrm{1}-E_{\mathrm{col}})a_{\mathrm{s}}{d}_{\mathrm{s}}^{b_{\mathrm{s}}+\mathrm{5}}{N}_{\mathrm{s}}^{\mathrm{2}}.\end{array}$$

To test the consistency of Eqs. (B9) and (B10), which were derived using different methods, we repeated the CNTRL and PHIL0.4 simulations but with Eq. (B9) replaced by the analytical solution for the change in number concentration from self-collection derived by Verlinde and Cotton (1993). The sensitivity of the results to this modification was found to be insignificant.

Graupel–graupel collisions are also parameterized in a similar manner. Since there is no graupel aggregation (collection efficiency of such collisions is assumed to be negligible), 100 % of the collisions are assumed to contribute to breakup:

$$\begin{array}{}\text{(B11)}& {\displaystyle \frac{\mathrm{d}{N}_{i_{\mathrm{gg}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1108}{a}_{\mathrm{g}}}{\mathrm{4}\times \mathrm{720}}}{\mathit{\pi }}^{\frac{\mathrm{1}-b_{\mathrm{g}}}{\mathrm{3}}}{\mathit{\rho }}^{\frac{\mathrm{2}+b_{\mathrm{g}}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{g}}^{\frac{-\mathrm{2}-b_{\mathrm{g}}}{\mathrm{3}}}{Q}_{\mathrm{g}}^{\frac{\mathrm{2}+b_{\mathrm{g}}}{\mathrm{3}}}{N}_{\mathrm{g}}^{\frac{\mathrm{4}-b_{\mathrm{g}}}{\mathrm{3}}},\text{(B12)}& {\displaystyle \frac{\mathrm{d}{Q}_{i_{\mathrm{gg}}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{836}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{48}\mathit{\rho }{\mathit{\rho }}_{\mathrm{g}}}}{a}_{\mathrm{g}}{d}_{\mathrm{g}}^{b_{\mathrm{g}}+\mathrm{5}}{N}_{\mathrm{g}}^{\mathrm{2}}.\end{array}$$

The value 1108 in Eq. (B11) is valid for b_s= 0.4 (Passarelli, 1978); in M05 b_s= 0.41 and b_g= 0.37; thus, adapting this value for both snow–snow (Eq. B9) and graupel–graupel (Eq. B11) collisions is a reasonable approximation.

Following the methodology of Sullivan et al. (2018b) in the TAKAH simulation, the number of fragments generated due to ice–ice particle collisions (F_BR) is

$$\begin{array}{}\text{(B13)}& F_{\mathrm{BR}}=\mathrm{280}(T-\mathrm{252}{)}^{\mathrm{1.2}}{e}^{-(T-\mathrm{252})/\mathrm{5}}.\end{array}$$

However, Takahashi et al. (1995) used 2 cm hail balls in their experiments; thus, to further include the influence of size in this formulation, we implement a size-scaled expression in the TAKAHsc simulation, assuming that F_BR depends linearly on D, decreasing to 0 at D= 0:

$$\begin{array}{}\text{(B14)}& F_{\mathrm{BR}}=\mathrm{280}(T-\mathrm{252}{)}^{\mathrm{1.2}}{e}^{-(T-\mathrm{252})/\mathrm{5}}{\displaystyle \frac{D}{D_o}},\end{array}$$

where D (in meters) is the size of the ice particle that undergoes fracturing and D_o= 0.02 m is the size of hail balls used by Takahashi et al. (1995).

The Phillips et al. (2017a) parameterization allows for the varying treatment of F_BR depending on the ice crystal type and habit.

$$\begin{array}{}\text{(B15)}& F_{\mathrm{BR}}=\mathit{\alpha }A\left(\mathrm{1}-\exp \left\{-{\left[{\displaystyle \frac{{\mathrm{CK}}_o}{\mathit{\alpha }A}}\right]}^{\mathit{\gamma }}\right\}\right),\end{array}$$

where $K_o=\frac{m_{\mathrm{1}}{m}_{\mathrm{2}}}{m_{\mathrm{1}}+m_{\mathrm{2}}}{\left(\mathrm{\Delta }{u}_{n_{\mathrm{12}}}\right)}^{\mathrm{2}}$ ,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\psi }=\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{3}},\\ {\displaystyle }& {\displaystyle }a=\mathit{\pi }{D}^{\mathrm{2}},\end{array}$$

where m₁ and m₂ are the masses of the colliding particles and Δu_n12 is the difference in their terminal velocities. The correction applied in Eq. (B7) is also adapted here to account for underestimates when u_n1≈u_n2. D (in meters) is the size of the smaller ice particle which undergoes fracturing, and α is its surface area. The parameterization was developed based on particles with diameters of 500 µm < D < 5 mm; however, Phillips et al. (2017a) suggest that it can be used for particle sizes outside the recommended range as long as the input variables to the scheme are set to the nearest limit of the range. C is the asperity-fragility coefficient, and ψ is a correction term for the effects of sublimation based on the field observations of Vardiman (1978). For cloud ice–snow, cloud ice–graupel, snow–graupel and snow–snow collisions,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }A=\mathrm{1.58}\times {\mathrm{10}}^{\mathrm{7}}\left(\mathrm{1}+\mathrm{100}{\mathrm{\Psi }}^{\mathrm{2}}\right)\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1.33}\times {\mathrm{10}}^{-\mathrm{4}}}{D^{\mathrm{1.5}}}}\right),\\ {\displaystyle }& {\displaystyle }\mathit{\gamma }=\mathrm{0.5}-\mathrm{0.25}\mathrm{\Psi },\\ {\displaystyle }& {\displaystyle }C=\mathrm{7.08}\times {\mathrm{10}}^{\mathrm{6}}\mathit{\psi }.\end{array}$$

The above parameters adapted from Phillips et al. (2017a) concern planar crystals or snow with a rimed fraction of Ψ < 0.5 that undergo fracturing: Ψ≤ 0.2 corresponds to lightly rimed particles, while Ψ≈ 0.4 represents highly rimed crystals/snow. The choice of the ice habit is based on particle images collected during the MAC flights which indicate the presence of needles and planar particles (O'Shea et al., 2017); needles are often considered secondary ice (Field et al., 2017). However, the M05 scheme does not explicitly consider habit and assumes spherical particles for all processes except sedimentation, for which the fall speed–diameter relationships are for non-spherical ice.

For graupel–graupel collisions, the parameters implemented in Eq. (B15) are somewhat different (Phillips et al., 2017a):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }A={\displaystyle \frac{a_o}{\mathrm{3}}}+\max \left({\displaystyle \frac{\mathrm{2}{a}_o}{\mathrm{3}}}-{\displaystyle \frac{a_o}{\mathrm{9}}}\left|T-\mathrm{258}\right|,\mathrm{0}\right),\\ {\displaystyle }& {\displaystyle }{a}_o=\mathrm{3.78}\times {\mathrm{10}}^{\mathrm{4}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{0.0079}}{D^{\mathrm{1.5}}}}\right),\\ {\displaystyle }& {\displaystyle }\mathit{\gamma }=\mathrm{0.3},\\ {\displaystyle }& {\displaystyle }C=\mathrm{6.3}\times {\mathrm{10}}^{\mathrm{6}}\mathit{\psi }.\end{array}$$

Finally, an upper limit for the number of fragments produced per collision is imposed and set to $F_{{\mathrm{BR}}_{\max }}=$ 100; this is the same for all collision types (Phillips et al., 2017a).

We estimate the production rate of fragments for cloud ice–snow collisions and cloud ice–graupel collisions using Eqs. (B2) or (B3) and one of the proposed formulations for F_BR above: $\frac{\mathrm{d}{N}_{i_{\mathrm{is}}}}{\mathrm{d}t}{F}_{\mathrm{BR}}$ and $\frac{\mathrm{d}{N}_{i_{\mathrm{ig}}}}{\mathrm{d}t}{F}_{\mathrm{BR}}$. For both of these collision types we assume that the cloud ice particles undergo breakup and the new smaller ice fragments remain within the same ice particle category. For snow–graupel collisions in which the snow particle is assumed to undergo fracturing, the production term $\frac{\mathrm{d}{N}_{i_{\mathrm{sg}}}}{\mathrm{d}t}{F}_{\mathrm{BR}}$ is added to the cloud ice category. In this case, mass transfer from the snow to the cloud ice category also occurs, but according to Phillips et al. (2017a), this is only 0.1 % of the snow mass that collides with graupel (B4). Snow–snow and graupel–graupel collisions are handled in the same way: $\frac{\mathrm{d}{N}_{i_{\mathrm{ss}}}}{\mathrm{d}t}{F}_{\mathrm{BR}}$ and $\frac{\mathrm{d}{N}_{i_{\mathrm{gg}}}}{\mathrm{d}t}{F}_{\mathrm{BR}}$ are added to the cloud ice number equation, while 0.1 % of $\frac{\mathrm{d}{Q}_{i_{\mathrm{ss}}}}{\mathrm{d}t}$ (Eq. B10) and $\frac{\mathrm{d}{Q}_{i_{\mathrm{gg}}}}{\mathrm{d}t}$ (Eq. B12) is added to the corresponding mass equation.