Simulations of primary and secondary ice production during an Arctic mixed-phase cloud case from the Ny-Ålesund Aerosol Cloud Experiment (NASCENT) campaign

The representation of Arctic clouds and their phase distributions, i.e., the amount of ice and supercooled water, influences predictions of future Arctic warming. Therefore, it is essential that cloud phase is correctly captured by models in order to accurately predict the future Arctic climate. Ice crystal formation in clouds happens through ice nucleation (primary ice production) and ice multiplication (secondary ice production). In common weather and climate models, rime splintering is the only secondary ice production process included. In addition, prescribed number concentrations of cloud condensation nuclei or cloud droplets and ice-nucleating particles are often overestimated in Arctic environments by standard model configurations. This can lead to a misrepresentation of the phase distribution and precipitation formation in Arctic mixed-phase clouds, with important implications for the Arctic surface energy budget. During the Ny-Ålesund Aerosol Cloud Experiment (NASCENT), a holographic probe mounted on a tethered balloon took in situ measurements of number and mass concentrations of ice crystals and cloud droplets in Svalbard, Norway, during fall 2019 and spring 2020. In this study, we choose one case study from this campaign that shows evidence of strong secondary ice production and use the Weather Research and Forecasting (WRF) model to simulate it at a high vertical and spatial resolution. We test the performance of different microphysical parametrizations and apply a new state-of-the-art secondary ice parametrization. We find that agreement with observations highly depends on the prescribed cloud condensation nuclei/cloud droplet and ice-nucleating particle concentrations and requires an enhancement of secondary ice production processes. Lowering mass mixing ratio thresholds for rime splintering inside the Morrison microphysics scheme is crucial to enable secondary ice production and thereby match observations for the right reasons. In our case, rime splintering is required to initiate collisional breakup. The simulated contribution from collisional breakup is larger than that from droplet shattering. Simulating ice production correctly for the right reasons is a prerequisite for reliable simulations of Arctic mixed-phase cloud responses to future temperature or aerosol perturbations.

1 Introduction

Given that the Arctic is the fastest-warming region on Earth, understanding the drivers of Arctic climate change, and in particular the role of clouds in this warming, has been of special interest (e.g., Serreze and Barry, 2011; Wendisch et al., 2017, 2019). Cloud phase is an important parameter in Arctic clouds as it influences both their radiative properties and their lifetime and thereby controls their climatic impact (e.g., Morrison et al., 2012). Mixed-phase clouds, i.e., clouds consisting of both ice crystals and supercooled liquid water droplets, can form at temperatures between approximately −38 and 0 °C. In the Arctic, these clouds frequently occur in every season and are especially common at low and middle levels (Shupe et al., 2006; de Boer et al., 2009; Shupe, 2011; Gierens et al., 2020). In order for ice to form in the mixed-phase cloud temperature regime, a special subset of aerosol known as ice-nucleating particles (INPs) is required. INPs act as a catalyst for water to freeze at temperatures above −38 °C (the homogeneous freezing temperature of water; Vali et al., 2015). Previous studies suggest that ice formation and thereby also the phase partitioning between liquid and ice is strongly affected by the availability of INPs (Jackson et al., 2012; Solomon et al., 2018; Norgren et al., 2018; Carlsen and David, 2022; Creamean et al., 2022). However, the degree to which the phase composition is determined by INPs is still an ongoing topic of research and investigation in field campaigns (e.g., Pasquier et al., 2022a).

There are different mechanisms through which INPs can nucleate ice, namely immersion, contact, condensation, or deposition freezing (Vali et al., 2015). Immersion freezing occurs when an INP is immersed in a liquid droplet and initiates freezing (Vali et al., 2015). Condensation freezing is similar to immersion freezing except that the INP acts as a cloud condensation nuclei (CCN) and nucleates ice as soon as bulk water condenses on its surface. Contact freezing occurs when an INP triggers a supercooled cloud droplet to freeze upon collision (Diehl et al., 2002; Yang et al., 2020). Lastly, deposition freezing occurs when water vapor directly nucleates into ice on an INP (Vali et al., 2015); however more recent studies have shown that deposition freezing may actually occur due to bulk water condensing in nanoscale cracks or pores, i.e., pore condensation and freezing (e.g., David et al., 2019b, 2020; Marcolli, 2014; Campbell and Christenson, 2018). Although condensation freezing is considered a separate ice-nucleation process, it is debatable whether it is microphysically different from deposition/pore condensation and freezing and immersion freezing (Kanji et al., 2017).

While these heterogeneous ice-nucleation processes, together with homogeneous nucleation, are commonly summarized as primary ice production processes, there has long been evidence for the existence of additional processes enhancing the number of ice crystals in clouds. Field measurements show that ice crystal number concentrations (ICNC) can significantly exceed INP concentrations (INPC) in clouds by up to 5 orders of magnitude (e.g., Auer et al., 1969; Beard, 1992; Cantrell and Heymsfield, 2005; Korolev and Leisner, 2020; Wieder et al., 2022a; Järvinen et al., 2022; Ladino et al., 2017). To explain this discrepancy, a number of so-called secondary ice production (SIP) processes for the multiplication of ice from existing frozen particles have been proposed. These processes include the collisional breakup of ice crystals (BR), rime splintering (RS; also called the Hallett–Mossop process), droplet shattering when freezing (DS), and sublimation fragmentation in subsaturated cloud regions (SF; e.g., Field et al., 2017; Korolev and Leisner, 2020). The existence of additional SIP processes has been proposed, but these have yet to be named and confirmed (Knight, 2012). RS occurs when liquid droplets (diameter < 13 µm) or drops (> 25 µm) rime on ice crystals but is only believed to be active in the temperature range between −8 and −3 °C (Hallett and Mossop, 1974; Mossop, 1978; Field et al., 2017). BR, caused by ice-ice collisions; DS, i.e., splinter production during the freezing of large droplets; and SF have been observed outside of this temperature range (Lauber et al., 2018; Korolev and Leisner, 2020; Keinert et al., 2020; Lauber et al., 2021, and references therein) and thus have the potential to significantly increase ICNC throughout the entire mixed-phase temperature range (Karalis et al., 2022).

Nevertheless, RS is generally the only SIP process that is part of the standard microphysics schemes available for both numerical weather prediction and climate models (e.g., Field et al., 2017; Zhao et al., 2021; Atlas et al., 2022). Many previous studies of model–observation comparisons therefore do not include additional SIP processes (Prenni et al., 2007; Luo et al., 2008; Barton et al., 2012; Young et al., 2017; Schemann and Ebell, 2020). These include, for example, the extensive modeling of Arctic stratocumulus clouds performed by Solomon et al. (2009) for cases in the Mixed-Phase Arctic Cloud Experiment (M-PACE) campaign that took place in 2004 (Verlinde et al., 2007). They tested several one- and two-moment schemes, all without SIP apart from RS, inside the Weather Research and Forecasting model (WRF; Skamarock et al., 2019) and found an underestimation of the number of small ice crystals. However, more recent studies have started including DS and BR to address the lack of SIP processes in models. Fu et al. (2019), for instance, studied a case from the same campaign using the Milbrandt and Yau microphysics scheme (Milbrandt and Yau, 2005) with different SIP configurations including tests for DS and BR. They found that BR did not contribute significantly in their case, while DS enhanced ICNC by roughly a factor of 2. Also, Sotiropoulou et al. (2021) and Georgakaki et al. (2022) implemented BR and DS, respectively, into the Morrison microphysics scheme (Morrison et al., 2009) in the WRF model following the parametrizations by Phillips et al. (2017, 2018). Sotiropoulou et al. (2021) found that including BR improved the representation of ICNC in the model considerably compared to field observations of Antarctic clouds. Georgakaki et al. (2022) found that including BR also increases the simulated ICNC in Alpine mixed-phase clouds and thereby leads to better agreement with observations, while DS did not significantly contribute to the simulated ICNC.

In our study, we use the WRF model including the modified Morrison microphysics scheme from Sotiropoulou et al. (2021) and Georgakaki et al. (2022) for a case study from a campaign in the Svalbard archipelago where information about the presence of SIP is available from field measurements (Pasquier et al., 2022b). The Ny-Ålesund Aerosol Cloud Experiment (NASCENT) study took place from September 2019 to August 2020 in Ny-Ålesund, Svalbard (78.9° N, 11.9° E; Pasquier et al., 2022a). One of its central objectives was to determine under which conditions INPs or secondary ice production (SIP) dominantly affects the phase partitioning inside Arctic mixed-phase clouds. For a deeper understanding of this dependency, we complement the observations performed with a regional modeling study here. The main objective of this study is to simulate ice production that is in agreement with observations, in terms of both numbers and processes. We are particularly interested in whether the model can simulate an ICNC as high as was observed when using measured INPC and added SIP processes.

A description of the observations used is given in Sect. 2, while modeling methods follow in Sect. 3. Section 4 comprises results and discussion before the conclusions are presented in Sect. 5.

2 Observations

We first present the relevant measurements from the campaign and then describe the meteorological situation of the chosen case study.

2.1 Instruments and measurements

In this study, we utilize the extensive suite of aerosol and cloud measurements conducted during the NASCENT campaign (Pasquier et al., 2022a), which took place over an entire year (fall 2019–2020) in and around the Norwegian research village Ny-Ålesund (78.9° N, 11.9° E). The measurements we use here are concentrations of aerosols suitable to act as CCN and INPs measured at the ground, as well as number and mass concentrations of ice crystals and cloud droplets measured in cloud from the HOlographic imager for Microscopic Objects (HOLIMO) on board the tethered balloon system HoloBalloon (Ramelli et al., 2020).

The CCN concentration (CCNC) was assessed from the concentration of aerosols larger than 70 nm, as particles of these sizes can be used as a proxy for CCN (Koike et al., 2019; Pasquier et al., 2022a).

To quantify the INP concentrations, aerosols were sampled through an inlet in a container installed at the balloon launch site using (1) the Horizontal Ice Nucleating Chamber (HINC; Lacher et al., 2017; Mahrt et al., 2018) and (2) a high-flow-rate liquid impinger (Bertin Coriolis-u), with subsequent analysis using the DRoplet Ice Nuclei Counter Zurich (DRINCZ; David et al., 2019a; Wieder et al., 2022b). Both techniques measure INPs in the immersion mode (Li et al., 2022). It is important to note that due to the efficiency of the liquid impinger, only aerosols larger than 500 nm could be assessed for their ice-nucleating ability in DRINCZ. However, previous studies have shown that aerosols larger than 500 nm make up the largest fraction of the INP population (DeMott et al., 2010; Mason et al., 2015). HoloBalloon itself, or more precisely the HOLIMO3B on board it, provided hydrometeor number and size concentrations measured at altitudes up to 850 m. The technical setup and pixel size restrict the detection of small particles to diameters larger than 6 µm and smaller than 2 mm (Ramelli et al., 2020). The holographic imager sampled a 15.5 cm³ volume of cloud at a rate of 6 Hz during the flights (Pasquier et al., 2022a, supplement). For this study, we use the hydrometeor number and size concentrations accumulated over 30, 60, or 300 s that were published in Pasquier et al. (2022c). All of the recorded hydrometeors larger than 25 µm were classified as either liquid droplets or ice particles using supervised machine learning (Touloupas et al., 2020). Meanwhile, all particles smaller than 25 µm were automatically classified as liquid droplets, as a reliable phase classification based on particle shape from HoloBalloon is limited to particles larger than this threshold (Lauber, 2020). Therefore, the retrieved ICNC is a lower estimate, strictly speaking. The phase-separated number size distributions were used to calculate the liquid water content (LWC) and ice water content (IWC). For the IWC, the effective mass–dimensional relationship reported in Cotton et al. (2013) was used. The associated uncertainties when using this approach are discussed in Heymsfield et al. (2010).

2.2 Case study

In this study, we focus on the 12 November 2019 cloud event from the NASCENT campaign, when there was observational evidence for a large contribution of secondary ice production to overall ice formation in the cloud (Pasquier et al., 2022a, b). Pasquier et al. (2022a, b) determined this by following the approach presented in Korolev et al. (2020), where the concentration of small pristine ice crystals with diameters < 100 µm is used as a way to identify regions where recent ice formation occurred. If the observed concentration of the newly formed ice crystals is larger than the INPC, it is likely that SIP was occurring. Still, this assumption includes the caveat that the actual INPC might be larger than measured since INPs in the contact mode were not assessed. For a more detailed description of the application during the NASCENT campaign, see Pasquier et al. (2022b).

On 12 November, a warm front influenced the weather around Ny-Ålesund (Fig. A1). The surface temperature varied between −3 and 0 °C, the dominant wind direction was southwesterly, and there was a persistent mixed-phase cloud lasting until around 21:00 UTC. As the cloud top rose from ca. 1300 to 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ during the day, the cloud top temperature decreased from around −11 to −14 °C. The cloud base varied between 200 and 600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, and there were several periods of precipitation resulting in about 2.4 mm of total precipitation (Pasquier et al., 2022a).

HoloBalloon performed three flights during 12 November 2019: 10:00–11:00, 12:15–14:00, and 14:45–17:00 UTC. The HoloBalloon measurement time spent at certain altitudes was not uniformly distributed; rather, the majority of the sampling time was spent close to the maximum altitude of the flight, leading to the highest robustness of the measurements inside the main body of the cloud. The estimated CCNC on 12 November of around 9 cm⁻³ is within a factor of 2 of the observed cloud droplet number concentration (CDNC) that varied between 5 and 15 cm⁻³, indicating that droplet formation was CCN limited (Motos et al., 2023). Such low concentrations were also commonly observed in earlier studies of CDNC in aerosol-limited regions (e.g., Mauritsen et al., 2011; Moore et al., 2013). The INP measurements made by DRINCZ during the day were fitted and reported as a function of temperature in Eq. (1) (Pasquier et al., 2022a, Fig. 7).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{n}_{\text{INP}}\left(T\right)\left[{\mathrm{m}}^{\mathrm{3}}\right]=& \mathrm{1000}\cdot \exp (-\mathrm{0.4146}\cdot (T\left[\mathrm{K}\right]-\mathrm{273.15})\\ & -\mathrm{12.4059})\end{array}\end{array}$$

The fit is several orders of magnitude lower than other parametrizations commonly used in weather and climate models that are often developed based on mineral dust or measurements from urban mid-latitude areas (e.g., Tobo et al., 2020), but it is on a similar order of magnitude to the INPC on other days during the NASCENT campaign (Li et al., 2022) as well as to results from field campaigns in other Arctic and remote coastal sites (e.g., Hartmann et al., 2020; McCluskey et al., 2018; DeMott et al., 2016; Sze et al., 2023).

The reported CCNC and INP fit are used to replace default values in the model. This permits better quantification of heterogeneous cloud particle formation, and also guides us towards better agreement of the modeled cloud particle concentrations with observations through the correct processes. Details about the model setup and implementation of measured aerosol properties follow in the next section.

Figure 1 The three nested domains D01, D02, and D03 in the WRF model are shown by the green, yellow, and red boxes, respectively. The associated horizontal resolutions are 15, 5, and 1 km.

3 Modeling setup and methods

All of the simulations presented here were conducted with the WRF model. To account for the different microphysics schemes investigated, two model versions were used, namely version 4.2.1 for simulations with the Milbrandt and Yau (MY) microphysics scheme (Milbrandt and Yau, 2005) and version 4.0.1 for simulations with the Morrison (Morr) microphysics scheme (Morrison et al., 2009; Skamarock et al., 2019). The simulations with the Morr scheme were conducted with this version of WRF so that the modified Morr scheme with the secondary ice production processes implemented by Sotiropoulou et al. (2021) and Georgakaki et al. (2022) could be used. We used a nested setup with three domains, where the outermost domain had a resolution of 15 km, the middle domain had a 5 km resolution, and the inner domain had a 1 km resolution. The geographical extent of the domains is shown in Fig. 1, where the innermost domain spans 100 km × 100 km. The number of vertical levels between the surface and 50 hPa was set to 172, 93 of which are below 3 km in altitude. This high vertical resolution was chosen to ensure that processes at cloud top and base were resolved in sufficient detail. All of the simulations were initiated on 11 November 2019 at 12:00 UTC and had a duration of 36 h with a time step of 30 s. The first 12 h of the simulations were used as a spinup. We initialized and nudged the model using reanalysis data from ERA5 on pressure levels at 00:00, 06:00, 12:00, and 18:00 UTC on the two outermost domains (Hersbach et al., 2018a, b). Grid nudging was performed for the two horizontal wind components, temperature, and specific humidity. No nudging was performed at the surface. As is by default in WRF, the nudging strength was 0.0003 s⁻¹, and nudging lasted for 60 min with a ramp-down at the end of the period. Longwave and shortwave radiation were treated by the CAM scheme (NCAR Community Atmosphere Model; Collins et al., 2004), and for boundary layer processes we used the YSU scheme (Yonsei University; Hong et al., 2006). The cumulus parametrization for deep and shallow convection on sub-grid scales was only turned on in the largest domain, and here we used the scale-aware Grell–Freitas ensemble scheme (Grell and Freitas, 2014).

The described settings for domain size, resolution, run duration, time step, nudging, radiation, and convection parametrizations are the same for all runs. Changes are only made inside the microphysics schemes and will be explained below.

3.1 Microphysics parametrizations

To investigate the sensitivity of cloud properties such as CDNC, ICNC, LWC, and IWC to different microphysical parametrizations and to the prescribed CCN and INPC, we perform similar sets of simulations with two different microphysics schemes that are described in more detail in the following subsections. Both schemes are double-moment schemes, but only the MY scheme predicts the CDNC from a prescribed CCNC and ambient supersaturation, while the Morr double-moment scheme uses a prescribed CDNC and is therefore actually only single-moment for cloud droplets. Additionally, MY has the largest number of hydrometeor classes, including cloud droplets, cloud ice, rain, snow, graupel, and hail, while Morr in its default setup does not include hail. The reasons for focusing most of the study on simulations with the Morr scheme were that the MY scheme failed to produce a suitable control simulation due to excessive graupel production when CCN and INP concentrations were adapted to observed values (see Sect. 4.3) and, more importantly, that we wanted to apply and test the new SIP implementation recently developed for the Morr scheme by Sotiropoulou et al. (2021) and Georgakaki et al. (2022).

Milbrandt and Yau scheme

In the MY scheme, CCNC and INPC are characterized individually. The default CCNC when there are negligible vertical motions or downdrafts is 200 cm⁻³ for continental aerosol and 80 cm⁻³ when maritime aerosol is selected. By default, continental aerosol is assumed, and the corresponding value is used across the entire domain. During updrafts, the number of activated CCN is parametrized following Cohard and Pinty (2000), which takes into account local temperature, updraft speed, and pressure.

The three possible pathways for primary ice production are heterogeneous ice nucleation, homogeneous freezing of cloud droplets, and freezing of rain drops (parametrized after Bigg, 1953). For heterogeneous ice nucleation, a switch enables either (1) the empirical parametrization from Meyers et al. (1992) for deposition, condensation freezing, and contact freezing parametrized following Young (1974) or (2) the parametrization from Cooper (1986) for deposition nucleation and no contact nucleation. As a default, option 1 (Meyers condensation/deposition nucleation and Young contact freezing) is active. Once primary ice has formed, the only SIP process included is RS (Hallett and Mossop, 1974). Ice mass may also increase from riming.

Cohard and Pinty (2000)Meyers et al. (1992)Young (1974)Bigg (1953)Pasquier et al. (2022a)Young (1974)

Table 1 Overview of simulations for 12 November 2019 with the Milbrandt and Yau microphysics (MY) scheme. The default scheme has two parametrizations for CCNC depending on the ambient vertical velocity w.

Meyers et al. (1992)Cooper (1986)Bigg (1953)Meyers et al. (1992)Cooper (1986)Bigg (1953)Meyers et al. (1992)Pasquier et al. (2022a)Meyers et al. (1992)Pasquier et al. (2022a)Meyers et al. (1992)Pasquier et al. (2022a)Meyers et al. (1992)Pasquier et al. (2022a)

Table 2 Overview of runs for 12 November 2019 with the Morrison double-moment microphysics scheme (Morr). The parameters changed during different model runs were CDNC, INPC in primary ice production processes, and lower threshold values for mixing ratios during rime splintering. In addition, we implemented the secondary ice processes in the Morr microphysics scheme introduced by Sotiropoulou et al. (BR; 2021) and Georgakaki et al. (DS; 2022).

Mod.thr. – modified thresholds.

To assess the effect of changing the given aerosol parametrizations towards a more realistic representation of the environmental conditions in our case, we compare two simulations performed with the MY scheme. The first simulation uses the default parametrizations only (MY_def). In the second simulation, both the CCNC and INPC are adapted to match the estimated values during 12 November (MY_adap) as follows: the CCNC is set to 9 cm⁻³ during updrafts, downdrafts and negligible vertical motions. Regarding heterogeneous ice nucleation, we chose to retain the default option (option 1) as it includes contact nucleation, which previous studies have shown is more important in Arctic mixed-phase clouds than deposition freezing (Morrison et al., 2005b, e.g.,) and is not active in option 2. We replace the formula for deposition/condensation freezing by Meyers et al. (1992) with the fit for INPC in the immersion mode given in Eq. (1) (Pasquier et al., 2022a). This is because immersion freezing is seen as the most common process for ice nucleation in mixed-phase clouds in the Arctic and in general in supercooled layers warmer than −27 °C (de Boer et al., 2010, 2011; Westbrook and Illingworth, 2011). For contact freezing, we use the parametrization by Young (1974), as no measurements of INPs in the contact freezing mode were conducted. For a further discussion of the treatment of contact nucleation, see the following section on the Morr scheme. Additionally, we turn off the freezing of rain drops in the model, as this process should only happen either because of an immersed INP or upon collisions with an ice particle (both of which are already accounted for). All important parameter differences between the two simulations are given in Table 1.

Morrison scheme

The Morr scheme is described in Morrison et al. (2009) and builds on an older version published in 2005 (Morrison et al., 2005a). It is a double-moment scheme for the hydrometeor species rain, ice, snow, and graupel but only single-moment for cloud droplets as the CDNC is a predefined number. It also contains a switch to include hail as a separate category, but this is not applied in this study.

In order for the Morr scheme to accurately represent the observed cloud properties for the correct physical reasons, several adaptations to the Morr scheme were tested and combined. The entire sequence of adaptations is shown in Fig. 2 and Table 2.

Figure 2 From default to final simulation: the different rows illustrate the steps from the default version of the Morrison scheme towards a version that represents cloud ice correctly. Adjusting the CDNC and primary ice parametrizations is done by matching concentrations to observed values of CCNC and INPC. The modifications regarding secondary ice processes are necessary in order to make SIP as efficient in the model as was observed in this case.

Sinclair et al. (2016)Young et al. (2019)Atlas et al. (2020, 2022)

Table 3 Mixing ratio thresholds in that have to be overcome (kg kg⁻¹) in order for rime splintering to happen in the Morrison microphysics scheme.

First, we adjusted the CDNC from the default value of 250 cm⁻³ (Morr₀) to 9 cm⁻³ (Morr₁), to match the observations within a factor of 2 (5–15 cm⁻³ observed; Pasquier et al., 2022a). Next, we adjusted primary ice production (Morr₂). Regarding processes involving INPs, the scheme contains the following: contact freezing parametrized after Meyers et al. (1992), deposition nucleation after Cooper (1986), and immersion freezing of cloud droplets and rain as parameterized by Bigg (1953). It is important to clarify here that Meyers et al. (1992) presented different parametrizations for deposition/condensation freezing and for contact freezing, and the formula used here is not the one mentioned during the description of the MY scheme above for condensation/deposition freezing. As the INP concentration during the NASCENT campaign was measured in the immersion mode, we use Eq. (1) retrieved in Ny-Ålesund (Pasquier et al., 2022a) to replace the formula by Bigg (1953) for the immersion freezing of cloud droplets. Along with this change, we introduce the following condition to limit ice nucleation by INPC: new ice crystals are only nucleated if supercooled cloud droplets are present and if the number of INPs at the given temperature exceeds the number of ice particles present. This limitation is necessary to prevent infinite ice nucleation, as the scheme is not aerosol aware, and INPs nucleated earlier in the cloud are not removed (Kärcher and Marcolli, 2021). Before, i.e., in the default scheme (Morr₀), the maximum number of ice crystals newly formed through immersion freezing was limited only by the number of available droplets. However, a limitation based on the number of ice crystals exists for deposition nucleation in the default scheme. As such, we effectively combine these two limitations (ice crystal and cloud droplet number) now for immersion freezing.

Based on the same rationale as for MY, we disable immersion freezing of rain and deposition nucleation but keep the default formula by Meyers et al. (1992) active for contact nucleation. Even though the extent to which contact freezing is actually occurring in the atmosphere is questionable (Ladino Moreno et al., 2013; Marcolli et al., 2016; Nagare et al., 2016), a sensitivity test with the Morr scheme (not shown) revealed that without contact freezing active and with immersion freezing parametrized after Pasquier et al. (2022a) as the only heterogeneous ice nucleation process, simulated ice crystal concentrations were far lower than the observed ones. Agreement could not be achieved by adding SIP the way it is done in this study. Therefore, we decided to keep contact freezing active (see also Sect. 4.5). In summary, Morr₂ and MY_adap contain comparable microphysical adjustments.

In Morr₃, we keep all of the previous adjustments from Morr₂ and add the SIP processes BR and DS following the parametrizations by Phillips et al. (2017, 2018) implemented into the Morr scheme for WRF by Sotiropoulou et al. (2021) and Georgakaki et al. (2022). As a result of implementing the BR parametrization following Phillips et al. (2017), the rimed fraction of snow and ice particles has to be prescribed, and here we use a rimed fraction of 0.4 in all simulations that include BR (Morr₃, Morr₄, and Morr₅). This corresponds to heavily rimed particles and was found to give reasonable results in Sotiropoulou et al. (2021) even though the value in nature is highly variable. We also note that in the scheme by Georgakaki et al. (2022), DS was allowed to be triggered by the freezing of rain drops in addition to during freezing induced by ice/snow/graupel–rain collisions (true SIP). Here, as we do not allow for immersion freezing of rain (but only cloud droplets), DS only includes events triggered by collisions.

Motivated by the results from Morr₃, we continued by adjusting the mixing ratio thresholds required for RS to occur (Morr₄) and finally the numbers of splinters created per RS event (Morr₅). We found that RS was not active in Morr₃ because the required mixing ratios for RS were never exceeded, and none of the newly implemented SIP mechanisms were triggered in a sufficient number to increase ICNC either (see Sect. 4.4, Appendix C). However, the required mixing ratios for BR and DS used by Sotiropoulou et al. (2021) and Georgakaki et al. (2022) were actually much lower than the ones given by Morrison et al. (2009) for RS. In that light, we adapted the thresholds for RS to occur. In particular, we lowered the snow and graupel water mixing ratio thresholds to match those required for BR during snow–snow and graupel–graupel collisions (10⁻⁸ kg kg⁻¹) and lowered the cloud liquid and rainwater mixing ratios to those required for DS induced by rain–ice crystal collisions (10⁻⁶ kg kg⁻¹). Earlier high-resolution modeling studies over the Southern Ocean and the Antarctic showed an improvement in the representation of low-level mixed-phase clouds in models when removing the RS thresholds in the Morr scheme (Sinclair et al., 2016; Young et al., 2019; Atlas et al., 2020, 2022; Sotiropoulou et al., 2021). Sinclair et al. (2016) and Young et al. (2019) removed only the thresholds for liquid mixing ratios, while Atlas et al. (2022) explicitly recommend completely removing all mass thresholds. The mixing ratio threshold values used in the different studies are given in Table 3.

A further discussion of the changes in RS efficiency follows in Sect. 4.5.

Figure 3 Simulated and radiosonde profiles of temperature (orange) and dew point temperature (blue) over Ny-Ålesund: here the radiosoundings at (a) 11:00 UTC, (b) 14:00 UTC, and (c) 17:00 UTC (solid lines) are shown along with the simulations MY_def (dotted) and Morr₀ (dashed) at the same times. During the majority of altitudes and times, the profiles for MY_def are not distinguishable, as its values are identical to those of Morr₀. For illustration of the meteorological situation, we additionally show the wind barbs from the sounding observations.

3.2 Model–observation comparison methodology

Before comparing observed and simulated cloud properties, we test the model's ability to simulate the observed environment by comparing meteorological variables. For that, ambient and dew point temperature from the three radiosondes launched on 12 November are compared with the simulations at the closest grid point and time step. To then compare the observed and simulated hydrometeor concentrations, we take the model grid point that is closest to the location of HoloBalloon and average the simulated cloud properties over the given flight. To ensure that our results were robust, we also calculated these averages including a few neighboring grid points, but this did not alter the results (not shown). To match the data acquisition from HoloBalloon, we only distinguish between liquid and ice and sum up cloud and rain droplets into one liquid category (Pasquier et al., 2022a). Hence, the total LWC and CDNC always refer hereafter to the combined liquid category. Similarly, all ice particles are merged into one ice category characterized by total IWC and ICNC; these include ice crystals, snow, and graupel in the Morr and MY schemes, as well as hail in MY. When comparing LWC, IWC, CDNC, and ICNC with the values measured in situ, we bin the observational data into 100 m intervals between 0 and 800 m altitude centered around 50 m, 150 m, etc. The cloud extended to higher altitudes than the tethered balloon was able to reach; therefore, a direct comparison is limited to the part of the cloud below 800 m. Finally, we compare the simulated and observed precipitation accumulated over the 24 h of 12 November 2019, as well as downward longwave radiation at the surface. Total precipitation was measured using a single-fenced gauge, while the downward longwave radiation measurements were performed and provided by the French–German AWIPEV Research Base (Maturilli, 2019). As solid- and mixed-phase precipitation is known to be underestimated when using single-fenced gauges (Wolff et al., 2015; Kochendorfer et al., 2017; Nitu et al., 2018), we adjust for this undercatch. In a simple estimate, we add 24 % to the measured total precipitation, which is the average undercatch reported for single-Alter-shielded gauges by Kochendorfer et al. (2017).

Figure 4 Total water content on 12 November as modeled by the final simulation, Morr₅. The macrophysical shape of the cloud agrees well with observations of radar reflectivity (see Pasquier et al., 2022b, Fig. 3c). The black lines show the observed cloud top as retrieved from radar measurements (Ebell et al., 2022), while the three flights performed by HoloBalloon are shown as grey lines.

Figure 5 (a) Total ICNC, (b) IWC, (c) CDNC, and (d) LWC on 12 November 2019 from the two simulations with the MY microphysics scheme. The different colors indicate the different balloon flights: flight 1 (10:00–11:00 UTC) in red, flight 2 (12:15–14:00 UTC) in violet, and flight 3 (14:45–17:00 UTC) in blue. The different line styles indicate different simulations (solid line for MY_def and dashed line for MY_adap). The bar diagrams show observations put together in 100 m bins between 0 and 800 m altitude, centered around 50, 150 m, etc. Note that the y axis is linear from the surface up to 800 m and logarithmic above to highlight the region where measurements are available. The vertical black line in (c) illustrates the constant CDNC used in the Morr simulations.

Figure 6 Same as Fig. 5 but for the simulations with the Morrison microphysics scheme and without showing droplet number concentration as this was predefined in the Morr scheme. Panels (a–c) show the simulations Morr₀–Morr₂, while (d–f) show the simulations Morr₃–Morr₅. The three balloon flights are marked by different colors; the simulations are distinguishable by line style. The profiles from Morr₂ (a–c, dotted line) and Morr₃ (d–f, solid line) look identical.

4 Results and discussion

Before comparing the influence of the various microphysical parameterization changes in the model, we verify the performance of the model nudging by comparing the simulated meteorological conditions with radiosonde observations (Fig. 3). It should be noted that the radiosonde observations are incorporated into the Global Telecommunication System and thus into ERA5 data, so no large differences should be expected. The simulated temperature profiles over Ny-Ålesund from the MY_def and Morr₀ simulations match the 11:00, 14:00, and 17:00 UTC radiosonde temperatures very well (mean deviation of 0.76 °C for both MY and Morr; see Fig. 3). The simulations also capture the dew point temperature profile well (mean deviation of 2.56 °C for MY and 2.55 °C for Morr), especially at lower altitudes up to 750 hPa where the cloud was present during the simulations (0.97 °C mean deviation for MY, 0.90 °C for Morr). With the simulated meteorological conditions validated when using the two microphysics schemes, we now assess the simulated macrophysical structure of the cloud before comparing the results from the MY and Morr schemes against each other. Finally, we discuss the impact of the changes within the respective microphysics schemes.

4.1 Cloud macrophysics

Before diving into the analysis of microphysical cloud properties, we show that the simulations and radar observations of the macroscopic evolution of the cloud on 12 November agree well (Fig. 4). Throughout the day, the observed and simulated cloud tops rise from ca. 1300 m at 06:00 UTC to around 2000 m at 18:00 UTC. This is true for the simulations with both microphysics schemes, and relative to the increase in cloud top height throughout the day, variations in cloud top height between the different simulations are small (Figs. 5 and 6). This indicates that the changes in cloud microphysical properties discussed in subsequent sections are primarily caused by differences in the specific microphysics schemes and not by accompanying changes in cloud structure.

The simulated IWC and ICNC extend to higher altitudes than the LWC (and CDNC) during all flights, indicating an ice-topped mixed-phase cloud (liquid only below ca. 1300 m vs. ice until 1450 m during the first flight and ca. 1800 m vs. 2000 m during the third flight; see Figs. 5 and 6). This result is in agreement with the Cloudnet classification (available from the Cloudnet data portal) that also identifies the cloud as ice topped (Ebell et al., 2022). Previous observational studies have frequently observed liquid-topped mixed-phase clouds in the Arctic (de Boer et al., 2009; Morrison et al., 2012, e.g.,). Therefore, the ice at cloud top in this case study may be a result of local effects (e.g., orography) or of the synoptic situation.

4.2 Microphysics scheme intercomparison

In order to assess the influence of the MY and Morr schemes on the simulated cloud microphysics, we compare the results of the simulations using the two schemes in their default configurations (MY_def and Morr₀, as described in Sect. 3.1). When compared with the HoloBalloon measurements, it becomes clear that the LWC is captured well by both simulations near the maximum flight altitude of HoloBalloon (see Figs. 5d and 6c). However, both MY_def and Morr₀ are unable to reproduce the LWC below 650 m, except for MY_def during flight 1 (Figs. 5d and 6c), while MY_def also underestimates CDNC in this altitude region (Fig. 5c). Simultaneously, MY_def and Morr₀ underestimate IWC, particularly during the first flight. While the IWC magnitude is wrong in the simulations, they both have an almost-constant value throughout the atmospheric layer considered, which is consistent with the observations. MY_def is in slightly better agreement with the observed ICNC than Morr₀ is, but both simulations maintain too few ice crystals near the surface. This comparison suggests that relative to the observations, the ice crystal aggregation may be too efficient in the default schemes, MY_def and Morr₀, as evidenced by the near-constant IWC and rapid decrease in ICNC towards the surface relative to the observations (Figs. 5a, b and 6a, b). The influence of constraining these two microphysical parametrizations with aerosol observations and the cause of these systematic biases in the default versions are discussed in the following sections.

4.3 Simulations using the Milbrandt and Yau microphysics scheme

As described in the previous section, MY_def struggles to maintain enough LWC and CDNC below 650 m (Fig. 5c and d). The underestimation of liquid at lower altitudes is not overcome when the CCNC is adapted to observed values in MY_adap, but the representation is improved slightly, especially during flight 1. Meanwhile, agreement between observed and simulated LWC is good around 750 m (Fig. 5d). Additionally, taking the lowest part of the sounding comparison below 900 hPa into account (Fig. 3), this indicates that at least during parts of the day the simulated water vapor pressure is too low just above the surface. Thus, the lifting condensation level, where dew point and absolute temperature are equal, is located too high in the model, which is visible from the soundings at 11:00 and 17:00 UTC (Fig. 3a and c). Nonetheless, at 14:00 UTC there is no such discrepancy between simulated and observed dew point at the lowest altitudes (Fig. 3b).

Compared to the CDNC and LWC, variations in the ICNC and IWC are much larger between the two simulations. In MY_adap, the ICNC increases from MY_def and even surpasses the observed ICNC (Fig. 5a). This is quite surprising, as the INPC is significantly reduced in MY_adap relative to MY_def. Upon further investigation, we find that the majority of the ICNC (and IWC) in MY_adap is comprised of graupel (see Fig. B1), which is similar to the simulation in Pasquier et al. (2022a) but was not present in the observations. To determine the major cause for this increase in graupel concentrations between MY_def and MY_adap, we systematically performed simulations with only certain changes included during each and found that the change in CCNC during updrafts is the determining factor for the increase in graupel amount. We therefore conclude that the CCNC changes lead to changes in the cloud droplet size distribution that make riming more efficient. In a simulation where all adaptions given in Table 1 are included except the CCNC during updrafts, which is kept parameterized following Cohard and Pinty (2000), the ICNC decreases and is lower than in MY_def, consistent with what we expected from lowering INP concentrations. In that case, the ICNC was 2 to 3 orders of magnitude lower than observed and not dominated by graupel (Fig. B2). These results show that MY_def simulated the observed ICNC reasonably well but for the wrong reasons, i.e., because of higher CCNC and INPC than were observed. Instead, we expect SIP to be the missing process needed to simulate the cloud microphysical properties correctly and for the right reasons. Before conducting an assessment of the effectiveness of SIP in the next section using the Morr scheme, we elaborate on the relevance of looking at the differences between flights as well as at agreement between simulated and observed precipitation and radiation.

Both the observations and simulations show variations in hydrometeor concentrations between flights representing different environmental conditions throughout the day. This increases the credibility of the study's results even though the scope is limited to one day. The inter-flight spread in MY_adap, however, reveals that this simulation does not correctly capture the changes in ICNC, IWC, and LWC between flight 1 and the subsequent flights.

Figure 7 (a) Accumulated precipitation on 12 November (24 h) divided into liquid (blue) and frozen (green) fractions and (b) downward longwave radiation at the surface (GLW) between 10:00 and 17:00 UTC. Values are shown for the two simulations with the MY scheme, MY_def and MY_adap; selected simulations with the Morr scheme; and observations (total precipitation shown in grey). As the observed precipitation was measured using a single-fenced gauge, there is likely an undercatch in the amount. Here, we adjust for that by adding 24 % following Kochendorfer et al. (2017, see Sect. 3.2). In the radiation panel (b), the orange lines illustrate the median, the boxes extend from the first to the third quartile of the data, and the whisker endpoints mark the farthest point inside 1.5 times the interquartile range from the box. The observed GLW was measured from the French–German AWIPEV Research Base (Maturilli, 2019).

In addition to hydrometeor profiles, we assess the precipitation that accumulated over 24 h. The lack of simulated CDNC near the surface and the stronger reduction in ICNC and IWC towards the surface compared to the observations lead to a strong underestimation of precipitation in both simulations (Fig. 7a).

Downward longwave radiation at the surface during the balloon flight times (10:00–17:00 UTC) is lower than observations in both simulations and lower in MY_adap than in MY_def (Fig. 7b). The overall underestimation of downward longwave radiation may be explained by a stronger simulated reduction in cloud water content towards the surface and thus a higher and colder cloud base on average, but the detailed explanation remains uncertain.

Overall, our results show that the simulations with the MY scheme struggle to represent observed precipitation and observed CDNC and LWC at lower altitudes. Even though simulated ICNC and IWC agree better with observations, this is due to a higher CCNC and INPC than observed in MY_def and due to graupel production simulated too high in MY_adap, whereas we expect a strong SIP to be the missing element.

4.4 Simulations using the Morrison microphysics scheme

In this section, we show the results from six simulations with the Morr scheme starting with the default version, Morr₀. Regarding ICNC, the maximum simulated value in Morr₀ is around 200 m⁻³ and thus, depending on flight number and altitude, at least 1 order of magnitude below the observations (Fig. 6a). Also, the liquid part of the cloud has a base at ca. 400 m, in contrast to the observations that measured liquid hydrometeors all the way to the surface (see Fig. 6c). This disagreement in LWC towards the surface is fixed by adapting the CDNC to the observed value of 9 cm⁻³ as was done in the Morr₁ simulation, which apart from the updated CDNC has the same settings as Morr₀. The reduced CDNC also results in an increase in the liquid precipitation simulated, from 0.12 mm in Morr₀ to 0.59 mm in Morr₁, during 24 h on 12 November (Fig. 7a). This can be explained by fewer and larger droplets that are more rapidly converted to rain and fall out. However, it is important to note that in both Morr₀ and Morr₁, the majority of the precipitation is simulated during the early morning hours of 12 November and not during daytime when the balloon flights took place (not shown). Meanwhile, as expected, the IWC and ICNC are not impacted and remain on the same order of magnitude as in Morr₀. Adapting the CDNC to the observed value results in a slight decrease in snowfall but not by as much as the rainfall increases, leading to an increase in total precipitation (Fig. 7a).

When adapting the primary ice production via the INP concentration fit retrieved from the campaign measurements (as described in Sect. 3.1) in Morr₂, the total ICNC and IWC decrease substantially, as expected. Ice (number and mass) no longer reaches the surface, and even above ca. 300 m the ICNC and IWC decrease by approximately 1 and 2 orders of magnitude, respectively. This results in the suppression of the simulated snowfall at the surface but also in a slight decrease in the total accumulated precipitation, to 0.48 mm (Fig. 7a).

Figure 8 Ice crystal number concentration tendencies of (a) contact freezing (CON; in blue) and immersion freezing (IMM; on top in red) and (b, c) collisional breakup (BR; in blue), rime splintering (RS; on top in red), and droplet shattering (DS; on top in grey). DS is shown partly transparent to visualize the fact that BR is generally active simultaneously. Panels (a) and (b) show the Morr₅ simulation, while panel (c) shows the SIP tendencies for the Morr₄ simulation.

We expect that the decrease in ICNC from modifying heterogeneous nucleation will be counteracted by increasing SIP in Morr₃ following Sotiropoulou et al. (2021) and Georgakaki et al. (2022). However, this does not reduce the large discrepancy between the simulated and observed ICNC, and the results of Morr₃ are in fact very similar to those of Morr₂ (Fig. 6a and d). To understand why implementing SIP in the simulation does not lead to substantial differences between Morr₂ and Morr₃, we examined the process tendencies from the model output in Morr₃ (similar to Fig. 8 where this is shown for Morr₄ and Morr₅). This allows us to quantify the contribution of the different processes, i.e., primary and secondary ice production processes, to the number of ice crystals formed.

In Morr₃ as well as in the following simulations, we see that immersion freezing is mainly active at cloud top where the coldest temperatures occur, while contact freezing dominates inside the cloud (shown for Morr₅ in Fig. 8a). Meanwhile, in terms of SIP, RS did not occur in Morr₃ over Ny-Ålesund, even though a large portion of the cloud was within the RS temperature range and riming was occurring. BR occurred occasionally in Morr₃ but produced fewer ice crystals than primary ice production, and DS occurred only sporadically but did not yield a substantial contribution to ice production (not shown).

Even though additional processes were implemented in Morr₃, they did not increase the ICNC compared to Morr₂, and the simulated ICNC still differed from the observations, as evidence of SIP was observed by HoloBalloon (Pasquier et al., 2022b). To investigate this discrepancy, we first adapted the required mixing ratio thresholds for RS to occur and then increased the number of splinters produced during RS in the following section (Sect. 4.5).

4.5 Role of secondary ice

Contrary to our expectations, activating the SIP processes added to the Morr scheme by Sotiropoulou et al. (2021) did not immediately increase the ICNC. To investigate this surprising finding, we first lowered the thresholds required for RS to occur. Specifically, we set the cloud liquid/rainwater and snow/graupel water mixing ratio thresholds to 10⁻⁸ and 10⁻⁶ kg kg⁻¹, respectively (see Table 3). The impact of lowering these thresholds is investigated in Morr₄ and results in a flight-dependent change in the ICNC. Because graupel concentrations were negligible compared to snow in simulations with the Morr scheme, we focus on splintering events involving snow. Then, for RS to be able to occur in the model at a given time and altitude, the temperature must be within the given range (−8 °C $<T<-\mathrm{3}$ °C), the snow water mixing ratio threshold must be overcome and, either the cloud liquid mixing ratio threshold or the rainwater mixing ratio threshold must be overcome. The altitude regions and flights in which the individual conditions and the joint condition are fulfilled are given in Table C1. This overview explains the differences in the simulated impact of the RS process during different flights well. We find that the simulated snow water mixing ratio and the temperature are most influential in defining the RS active region, whereas the liquid water mixing ratio thresholds do not substantially limit the region further. As a result, during flights 1 and 2, the simulated ICNC increases by a factor of up to 2 and up to 10, respectively, while during flight 3 hardly any change in ICNC is simulated (see Fig. 6d, simulations Morr₃ and Morr₄). The latter is because the simulated snow water mixing ratio does not overcome the threshold at any altitude during the third flight. It should be noted that this is in contradiction to the observations, as HOLIMO measured the most SIP during flight 3. Along with RS, BR also became active during flights 1 and 2. Meanwhile, Morr₄ only showed a very minor contribution to the increase in ICNC from DS, and only during flight 2 (Fig. 8c). We discuss the relative importance of the different SIP mechanisms after presenting the results of the last adaption made in the scheme regarding the number of splinters produced during RS.

Previous studies have shown that the efficiency of the RS process increases in the presence of large cloud droplets (e.g., Mossop, 1978). As relatively large cloud droplets were observed during this case (Pasquier et al., 2022b, a), here we assess the impact of increasing the number of splinters produced by the RS process by a factor of 20 in the Morr₅ simulation (see Sect. 3.1). We justify this step by considering that cloud droplets in Arctic pristine conditions may be larger and may lead to more splinters.

Even though the mixing ratio thresholds are not modified between Morr₄ and Morr₅, increasing the number of splinters generated influences the total ICNC during all flights. This increase in ICNC leads to a subsequent increase in BR, which in combination with the initial increase in ICNC from RS leads to a total increase in ICNC by up to 3 orders of magnitude, with BR making up a large share of that increase (Fig. 8b). This effectively leads to a “cascading” process of ice production, as was observed in field observations (Lawson et al., 2015; Pasquier et al., 2022b).

When comparing the observed and simulated SIP mechanisms, the observations indicate a more prominent contribution of SIP from DS (Pasquier et al., 2022a). In contrast, Morr₅ simulates that DS is not a determining factor in the increase in ICNC and is only active at altitudes below 500 m/near cloud base (Fig. 8b). However, as the in situ observations were made near cloud base, the fact that Morr₅ also simulates a contribution from DS in this region indicates that even though the model underestimates the contribution from DS, it is simulating it in the correct place. Additionally, it should be noted that near cloud base, where DS was observed and simulated, the temperatures were within the RS range. This, in combination with the observed frozen/fragmented droplets and rimed columns/lollipop ice (Pasquier et al., 2022a), makes it difficult to disentangle the true contributions of RS and DS from the measurements.

Nevertheless, our result which showed that DS had a minor effect on the simulated ICNC is in agreement with the modeling study of wintertime alpine mixed-phase clouds by Georgakaki et al. (2022), which was conducted with the same modified Morrison microphysics scheme in WRF. This suggests a bias toward low DS activity in this scheme. As a possible key to alter this bias, two studies using the Morr microphysics scheme in global models suggest that if a different size representation (i.e., a bin representation for the radius instead of a bulk representation of the particle mass) is used for cloud droplets during the treatment of DS, the magnitude of DS increases (Zhao et al., 2021; Sotiropoulou et al., 2024).

BR played the largest role in increasing ICNC, not only in our case but also in those of Georgakaki et al. (2022) and Sotiropoulou et al. (2021). However, Georgakaki et al. (2022) found that in their case ice crystals falling from above were important in initiating BR, which was not the case here. Sotiropoulou et al. (2021) also found that even with RS completely deactivated, BR alone was able to represent the observed ICNCs in Antarctic summer clouds. At the same time, they admit that primary ice production might have been overestimated in their setup.

In our simulations with constrained CCNC and INPC and thus a more realistic representation of primary ice production, we chose to reduce the thresholds for RS instead of removing them, as Sotiropoulou et al. (2021) hypothesized that removing all mass thresholds for RS overestimated RS. As previously stated, when the number of splinters is also increased relative to the default parametrization, our simulation results match the observed ICNC very well (Fig. 6d). This result is in accordance with Young et al. (2019), who also found the best agreement between model and observations when removing the RS liquid mixing ratio threshold and when making RS 10 times more efficient in Arctic clouds. This may be due to the high concentration of large droplets observed in this case and typically found in the Arctic. When comparing the patterns in spatial and temporal extent of the SIP active regions for BR, RS, and DS (Fig. 8b and c), it becomes clear that the higher number of splinters from RS in Morr₅ also strongly affects the activity of BR and DS. For example, during flight 3 where RS was not simulated in Morr₄ (Fig. 8c), SIP still became active enough to represent the observed ICNC in Morr₅ after increasing the number of splinters produced (Fig. 6f, dotted blue line). Overall, the activity of RS, BR, and DS in time changes, e.g., due to the advection of ice particles between grid boxes, and a higher ICNC is maintained, again favoring more SIP in turn. Finally, our results show that in situations with constrained CCNC and INPC, maintaining some thresholds for RS works well when the number of splinters is also increased relative to the default parametrization.

As for the MY simulations, we assess the accumulated precipitation over 24 h, and here, we find a clear improvement in the representation in Morr₅, especially of solid precipitation (Fig. 7a). While highly underestimated in Morr₀, now total precipitation is even overestimated, likely due to an amount of supercooled drizzle that is too high. However, the phase assessment of precipitation from observations is uncertain and of a qualitative type only, as it is based on the notes of the campaign crew, which do not cover the whole 24 h period. Looking at the downward longwave radiation during the flights, the variations in the median are small between the simulations Morr₁, Morr₂, and Morr₅, which all show lower values than Morr₀ (Fig. 7b). This hints at the fact that the CDNC modification is influencing the radiation the most in our case. However, the reason for the underestimation compared to observations remains unclear.

Overall, Morr₅ clearly represents the best match with observations of ICNC and IWC, even though the relative importance of BR and DS does not necessarily reflect the observations. It should also be noted that contact freezing is still active in this simulation. Assuming that this process is not as important in nature as was simulated, the strength of SIP would need to be increased further. Morr₅ is the only simulation that produces as high of an ICNC as was observed, and it still represents LWC well and better than the default Morr scheme (Morr₀). In addition, it performs best in simulating enough precipitation. The combination of implemented and enhanced SIP processes is able to reproduce the observed ICNC and IWC during all flights, although the environmental conditions and the amount of IWC produced through primary ice production vary. Therefore, we consider Morr₅ a successful attempt to represent an Arctic mixed-phase cloud with a realistic distribution of primary and secondary ice production.

Figure 9 Temperature profiles from the model simulations averaged over the three flight times on 12 November. Here, the profiles from simulation Morr₄ are shown, but the choice of simulation does not play a role as the profiles do not vary across simulations. Temperatures below −38 °C, i.e., suitable for homogeneous freezing, are reached above ca. 6300 m (not shown). The vertical dashed lines and shading indicate the temperature range where rime splintering is permitted (−8 to −3 °C), while the colored horizontal dashed lines indicate the maximum flight altitudes.

5 Conclusions

This study shows that generalized out-of-the-box cloud microphysics schemes, i.e., MY and Morr, fail to correctly represent the vertical structure of ice and liquid water content of Arctic mixed-phase clouds. While these schemes do reproduce the observed maximum values reasonably well, we find that this occurs for the wrong reasons due to compensating errors. In the default schemes, the number of aerosols available for nucleating cloud particles is often unrealistically high for Arctic conditions, and the overly efficient primary ice production compensates for the lack of secondary ice production. As shown here, in the pristine Arctic environment of Ny-Ålesund, Svalbard, the low-level mixed-phase cloud observed on 12 November 2019 is well represented in the WRF model by the default MY and Morr schemes in terms of macrophysical shape and water content, but the ICNC and the amount of precipitation from the cloud are underestimated. When the concentrations of aerosols suitable to nucleate cloud particles are adapted to observed values in the MY scheme, the ICNC increases but is dominated by graupel, which was not observed but was already pointed out in Pasquier et al. (2022a). This increase in graupel is caused by the lowered CCNC during updrafts, indicating that the MY scheme fails to represent ice production accurately in environments with very low CCNC. Meanwhile, the underestimation of liquid water towards the surface and of precipitation is not altered when aerosol concentrations are changed. In the Morr scheme, however, the simulation matches the observed ICNC and precipitation well: when the concentrations of cloud nucleating aerosols are adapted to observed values, rime splintering is made more efficient, and further secondary ice processes are added. As the simulated primary ice production includes contact freezing, which is not observationally constrained and is generally thought to be of limited importance, it must be considered an upper bound for the plausible primary ICNC contribution in this case. If we instead assume that contact freezing was significantly less important in our case than was assumed in the simulations, or even not active at all, the already greatly enhanced secondary ice production would have to be even further enhanced in order for agreement between simulated and observed ICNC to be achieved.

The lack of modeled secondary ice production in the default Morr scheme is due to both missing processes (in our case mainly collisional breakup) and snow, cloud liquid, and rainwater mixing ratio thresholds that are too high for rime splintering to occur. Based on the earlier studies by Atlas et al. (2020, 2022); Sotiropoulou et al. (2021, 2020) and the finding in our case that the snow, cloud liquid, and rainwater mixing ratios never reached the thresholds required for rime splintering to happen in the default scheme, we recommend lowering these mixing ratio thresholds. The fact that the occurrence of rime splintering was needed in order to activate further SIP processes and ultimately represent observed ICNC further supports this recommendation. Otherwise, the thresholds restrict the occurrence of rime splintering to clouds with a higher ice and liquid water content than is often observed in the Arctic. Additionally, we increase the number of splinters produced and find that an increase by a factor of 20 yields good agreement with observed ice crystal number concentrations. Based on this and on the findings of Young et al. (2019), we hypothesize that this is justified in clean Arctic environments with relatively large droplets, but further laboratory and field studies relating the number of splinters produced by rime splintering to both mixing ratios and droplet sizes are required to confirm or reject this hypothesis. In addition to the ICNC profiles, the simulated precipitation is also considerably improved when aerosols and SIP processes are represented more realistically. The adapted schemes show a lower downward longwave radiation than default schemes, and the changes are likely mainly related to CCNC/CDNC changes. Thus, the misrepresentation of ice production might not lead to a bias in how much heat is trapped between the cloud and the surface, but due to an overall underestimation compared to observations, conclusions regarding downward longwave radiation remain uncertain.

Even though the changes made between the default simulation Morr₀ and the final simulation Morr₅ have been validated using the case study on 12 November with observational evidence for strong SIP, the applied modifications are meant to be of a general nature and should not hinder the application to cases where no SIP was observed. Also, the performance across different microphysical conditions between the flights already suggests that the modified scheme is applicable in different conditions. The correct representation of hydrometeor profiles and precipitation for the right reason is a prerequisite for model simulations to be used as baselines to investigate cloud responses to aerosol perturbations and future warming.

Appendix A: Weather map

The following figure taken from Pasquier et al. (2022a, Fig. S2) illustrates the meteorological situation around Svalbard on the morning of 12 November 2019.

Figure A1 Map showing the synoptic situation with the estimated location of the warm front on 12 November at 06:00 UTC (figure taken from Pasquier et al., 2022a, Fig. S2, produced using MetCoOp Ensemble Prediction System (MEPS) model data from Hellmuth and Hofer, 2019). The relative humidity and the temperature at 1000 hPa are indicated by green shading and colored lines, respectively. The red triangle shows the location of Ny-Ålesund.

Appendix B: Impact of adapting updraft CCNC in the MY scheme

As mentioned in Sect. 4.3, the surprisingly high amount of graupel in MY_adap (see Fig. B1) was not produced when all changes between MY_def and MY_adap were made, with the exception of fixing the CCNC during updrafts to the observed value. Instead, the default parametrization by Cohard and Pinty (2000) was used here. The results of this simulation are marked as MY_CP00 in Fig. B2.

Figure B1 Number concentrations of individual hydrometeor categories in simulation MY_adap during the third flight.

Figure B2 ICNC for the simulations MY_def, MY_adap, and MY_CP00, which is identical to MY_adap apart from the CCNC during updrafts. The ICNC is much lower when CCNC during updrafts is parametrized according to Cohard and Pinty (2000) than when it is set to 9 cm⁻³.

Appendix C: Impact of RS mixing ratio threshold changes on simulated ICNC

Whether lowering the required mixing ratio thresholds from simulation Morr₃ to Morr₄ has an effect on the simulated ICNC during one flight depends on whether the new thresholds are overcome during the time of this flight. To reveal the regions where the lowered mixing ratio thresholds for RS in Morr₄ actually enable RS, we show the altitude regions where the individual thresholds are overcome and where the joint condition (see Sect. 4.5) is met in Table C1. During flight 1, the flight-time-averaged snow mixing ratio in Morr₄ just overcame the updated threshold (by less than 50 %) at altitudes between 350 and 620 m. At the same time, the temperature criterion was met between 100 and 700 m in altitude, and the rain mixing ratio exceeded the threshold above 370 m. Thus, the altitude region where RS involving snow and rain was active during flight 1 is located between 370 and 620 m. RS involving snow and cloud droplets did not occur during flight 1 as the cloud liquid mixing ratio threshold was not overcome in the altitude region where the amount of snow was sufficiently high. During flight 2, both the snow and rain mixing ratios clearly exceeded the thresholds at all altitudes. Therefore, the impact of lowering the required mixing ratio thresholds is mostly captured during the second flight. The ice part of the cloud now extends to the surface as observed. Similar to during flight 1, the temperature profile leads to RS primarily increasing the ICNC towards cloud base where the cloud is warm enough for RS to be active (Fig. 9). During the third flight, however, the flight-averaged snow water mixing ratio was even lower than in the first flight and did not reach the threshold. This explains the lack of difference in ICNC between Morr₃ and Morr₄ for the third flight.

Table C1 Regions where criteria for RS are met in Morr₄. In order for RS to happen, the threshold for the snow mixing ratio q_snow has to be overcome, the temperature must be inside the given range, and the cloud water mixing ratio q_{cloud  liquid} and/or rainwater mixing ratio q_rain must overcome the threshold. The combination of these limitations results in the RS active region given in the lowest row.