RaFSIP: Parameterizing Ice Multiplication in Models Using a Machine Learning Approach

Accurately representing mixed‐phase clouds (MPCs) in global climate models (GCMs) is critical for capturing climate sensitivity and Arctic amplification. Secondary ice production (SIP), can significantly increase ice crystal number concentration (ICNC) in MPCs, affecting cloud properties and processes. Here, we introduce a machine‐learning (ML) approach, called Random Forest SIP (RaFSIP), to parameterize SIP in stratiform MPCs. RaFSIP is trained on 16 grid points with 10‐km horizontal spacing derived from a 2‐year simulation with the Weather Research and Forecasting (WRF) model, including explicit SIP microphysics. Designed for a temperature range of 0 to −25°C, RaFSIP simplifies the description of rime splintering, ice‐ice collisional break‐up, and droplet‐shattering using only a limited set of inputs. RaFSIP was evaluated offline before being integrated into WRF, demonstrating its stable online performance in a 1‐year simulation keeping the same model setup as during training. Even when coupled with the 50‐km grid spacing domain of WRF, RaFSIP reproduces ICNC predictions within a factor of 3 when compared to simulations with explicit SIP microphysics. The coupled WRF‐RaFSIP scheme replicates regions of enhanced SIP and accurately maps ICNCs and liquid water content, particularly at temperatures above −10°C. Uncertainties in RaFSIP minimally impact surface cloud radiative forcing in the Arctic, resulting in radiative biases under 3 Wm−2 compared to simulations with detailed microphysics. Although the performance of RaFSIP in convective clouds remains untested, its adaptable nature allows for data set augmentation to address this aspect. This framework opens possibilities for GCM simplification and process description through physics‐guided ML algorithms.


Introduction
In the big data and machine learning era, global climate models (GCMs) remain an indispensable tool for predicting how the Earth system will respond to rising greenhouse gas concentrations (Balaji et al., 2022;Irrgang et al., 2021).GCMs typically represent atmospheric processes with a horizontal grid spacing of about 50-100 km.However, at this coarse resolution, small-scale processes, such as those associated with clouds, cannot often be explicitly resolved.These complex subgrid scale processes interact with the resolved model scales only through parameterizations tuned in a way to improve the agreement with theory or observations (Hourdin et al., 2017).Parameterization schemes are a necessary but non-trivial part of climate modeling, introducing uncertainties in future climate projections.These uncertainties frequently result in persistent model biases, especially over • A random-forest parameterization for secondary ice production is developed using outputs from a 10-km horizontal grid spacing simulation • Cloud phase partitioning agrees within a factor of 3, with radiative biases below 3 Wm 2 compared to the detailed microphysics simulation • The scheme can be adjusted to coarser resolutions typical of climate models without losing computational efficiency and numerical stability Supporting Information: Supporting Information may be found in the online version of this article.
Machine learning (ML)-based parameterizations are a promising and computationally efficient approach increasingly used in climate science to replace, supplement or speed up conventional parameterizations (e.g., Brenowitz & Bretherton, 2019;Gentine et al., 2018;Grundner et al., 2022;Han et al., 2020;Mooers et al., 2021).ML parameterizations are sometimes developed based on high-resolution three-dimensional (3D) simulations, in which subgrid processes are either explicitly resolved (Brenowitz & Bretherton, 2018;Yuval et al., 2021;Yuval & O'Gorman, 2020) or parameterized using conventional schemes (Grundner et al., 2022;O'Gorman & Dwyer, 2018).In this approach, high-resolution data are often coarse-grained to match the lower spatial resolution of the GCM grid (e.g., Brenowitz & Bretherton, 2019).Advances in both computational tools and data assimilation methods have now allowed for global observations to be integrated into ML training data sets (Schneider et al., 2017), yet limitations concerning the limited spatiotemporal coverage of observational data still need to be addressed (Irrgang et al., 2021).
Neural networks (NNs) and random forests (RFs) are two ML methods that have been used for developing parameterizations for GCMs.RF-an ensemble learning algorithm consisting of multiple decision trees-ensures the preservation of physical properties by taking the average predictions of subsets of the training samples (Breiman, 2001).This conservative behavior leads to more stable simulations when coupled with GCMs compared to NNs (Han et al., 2020;O'Gorman & Dwyer, 2018), which can deviate from the training data pool.Notably, RF parameterizations have demonstrated successful emulation of conventional convective schemes in GCMs (O'Gorman & Dwyer, 2018).Yuval and O'Gorman (2020) applied RF parameterization using highresolution outputs from an idealized simulation that were coarse-grained, achieving effective replication of climate patterns across different grid spacings.In contrast, NN parameterizations, although requiring less memory in GCMs, can introduce numerical instabilities and result in model runaway or climate drift (Brenowitz & Bretherton, 2018, 2019;Rasp et al., 2018).To ensure stable long-term simulations, post-prediction adjustments and accurate processing of the training data set are necessary (Rasp, 2020;Yuval et al., 2021).While offline comparisons favored NN parameterizations, both RF and NN algorithms satisfied physical constraints (energy conservation and non-negative surface precipitation) and demonstrated comparable performance when coupled with the dynamical core of a coarse-resolution aqua-planet GCM (Yuval et al., 2021).Additionally, RFs and NNs showed similar performance in predicting the probability density function of surface solar irradiance under 3D cumulus clouds (Gristey et al., 2020).
Parameterizations of cloud microphysical processes and aerosol-cloud interactions are perhaps one of the most challenging aspects of climate simulations (e.g., Liu et al., 2023;Morrison et al., 2020).They are also computationally expensive because they require a large amount of information to be carried around in simulations (e.g., cloud microphysical and aerosol quantities and their interactions) and the parameterizations themselves also may require considerable computational effort.A few studies to date have tried to simplify cloud microphysical processes using ML-based approaches, focusing on warm cloud microphysical processes.For instance, Seifert and Rasp (2020) employed a Monte Carlo superdroplet simulation to train neural network (NN) models for parameterizing autoconversion, accretion, and self-collection rates in two-moment schemes.Chiu et al. (2021) developed NN-based parameterizations for autoconversion and accretion rates, incorporating in-situ observations of droplet size distributions and employing the stochastic collection equation with bin microphysics to separately account for cloud and drizzle water contents.Although these ML approaches were evaluated offline, Gettelman et al. (2021) successfully replaced the computationally expensive bin microphysical scheme used in a GCM with ML emulators to predict autoconversion and accretion tendencies, achieving comparable accuracy with reduced computational costs.However, there remains a research gap regarding ML parameterizations specifically tailored to represent ice cloud microphysics.
This study focuses on developing an RF-based approach for parameterizing secondary ice production (SIP) in stratiform mixed-phase clouds (MPCs) which are prevalent in polar regions (de Boer et al., 2009;Shupe et al., 2006), sustained by a complex interplay of microphysics, dynamics, radiatively driven turbulence, as well as surface heat and moisture fluxes (Morrison et al., 2012).Misrepresentation of mid-and high-latitude MPCs in GCMs has been shown to increase the spread of predicted cloud feedbacks in the recent Climate Model Intercomparison Project (Phase 6-CMIP6; Murray et al., 2021;Zelinka et al., 2020).Furthermore, the ice and liquid phase partitioning in GCMs can significantly impact Arctic amplification, highlighting the necessity of revisiting the microphysical parameterizations associated with Arctic MPCs (Tan & Storelvmo, 2019).
Ice multiplication, also known as SIP, can be an important source of ice particles in MPCs, as it can rapidly increase the pre-existing ice crystal number concentrations (ICNCs) through a number of collisional processes (Field et al., 2017;Korolev & Leisner, 2020).The most frequently acknowledged SIP mechanisms are the Hallett-Mossop rime splintering (HM; Choularton et al., 1980;Hallett & Mossop, 1974), ice-ice collisional break-up (BR; Takahashi et al., 1995;Vardiman, 1978) and droplet-shattering (DS; Choularton et al., 1980;James et al., 2021;Kleinheins et al., 2021).Observations of MPCs in the Arctic region have shown high ICNCs that greatly surpass ice nucleating particle (INP) concentrations (Wex et al., 2019), particularly at temperatures above 25°C (Luke et al., 2021;Pasquier et al., 2022;Rangno & Hobbs, 2001;Schwarzenboeck et al., 2009).This suggests that SIP might be a prevalent process in the moderately cold polar conditions.Along with observations, modeling studies across scales have also emphasized the significance of incorporating SIP parameterizations in regional and global climate model simulations for an accurate representation of the phase partitioning and radiative properties of Arctic MPCs (e.g., Fridlind & Ackerman, 2018;Sotiropoulou et al., 2020Sotiropoulou et al., , 2021;;Zhao et al., 2021;Zhao & Liu, 2021, 2022).Recent physically based SIP parameterizations tested in GCMs were developed by pooling experimental observations and further considering the physics of collisions or applying correction factors to account for the simplified laboratory setups (Phillips et al., 2017(Phillips et al., , 2018)).Zhao et al. (2021) and Zhao and Liu (2022) developed a hybrid-bin framework for improved representation of SIP, but this approach can be computationally expensive (Sotiropoulou et al., 2024).These advanced SIP formulations require that the host model explicitly describes a number of ice (cloud ice, snow and graupel) and liquid species (cloud droplets and raindrops) as well as interactions among them (i.e., aggregation, collection and riming).This complexity poses challenges for implementing inline prognostic SIP parameterizations in largescale models lacking detailed microphysics modules.Additionally, the inherent nonlinearity of processes like SIP, if explicitly implemented at coarse resolutions, can introduce significant biases.This, in turn, is a major contributing factor to models underestimating ICNCs by orders of magnitude in the presence of SIP.
Seeking to address the challenges posed by implementing complex microphysics at coarse spatial scales, we propose a new approach toward parameterizing SIP in stratiform MPCs, which we call the RaFSIP (Random Forest SIP) parameterization.RaFSIP is an ML data-driven parameterization trained on comprehensive model outputs from mesoscale model simulations that account for all three SIP processes: HM, BR and DS.Using only a few thermodynamic and microphysical input variables extensively described in models, RaFSIP expresses the effect of SIP either directly by predicting the corresponding SIP rates or through the use of the ice-enhancement factor (IEF), which is a multiplication factor applied to primary ice crystals.The ML approach presented in this paper demonstrates a way toward model simplification that maintains important aspects of higher resolution simulations, but remains compatible with simplified model frameworks.

Data and Methods
The reference simulation used to develop the RF framework is derived from the mesoscale Weather Research and Forecasting (WRF) model.The WRF outputs from the target simulation were used to train two different versions of the RaFSIP parameterization: the first version predicts the effect of SIP through the use of a multiplication factor applied to the primary ice production rates (hereafter denoted as RaFSIPv1), while the second version is trained to predict directly the production rates of secondarily formed ice particles (hereafter referred to as RaFSIPv2).The workflow followed to develop each parameterization is summarized in Figure 1.

WRF Regional Climate Simulations
The target simulation was generated using WRF version 4.0.1,including augmented cloud microphysics that considers the effects of BR and DS mechanisms introduced by Sotiropoulou, Vignon, et al. (2021) and Georgakaki et al. (2022), in addition to the default HM parameterization.A regional climate configuration was employed (Abdelwares et al., 2018), where a 50-km horizontal grid spacing parent domain covering the Arctic region contains a single-way nested domain of 10-km horizontal grid spacing.Figure 1 shows the map of the two domains centered over Ny-Ålesund using a polar-stereographic projection, that is suitable for high-latitude WRF domains.The outermost domain consists of 148 × 148 grid points, while the innermost one contains 301 × 301 grids.In the hybrid terrain-following mass coordinate system adopted by WRF, we used 40 vertical levels up to a model top of 50 hPa (i.e., ∼20 km), with the standard grid spacing generated automatically by WRF.
The initial and lateral boundary conditions are taken from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyzes data set (ERA5; Hersbach et al., 2020), with a 31-km horizontal grid spacing.The lateral forcing at the edge of the parent domain was updated every 6 hr.To ensure realistic synoptic dynamics in the model the outer domain of WRF has been nudged toward ERA5 reanalysis for zonal and meridional wind speed, with a relaxation time scale of 6 hr.Static fields come from default WRF pre-processing system data sets with a resolution of 30'' for both the topography and land use fields.The WRF physics options include the Rapid Radiative Transfer Model for General Circulation Models (RRTMG) radiation scheme to parameterize both shortwave and longwave spectra, the local Mellor-Yamada-Nakanishi-Niino Level 2.5 (MYNN; Nakanishi & Niino, 2006) scheme with its associated surface layer scheme for the planetary boundary layer (PBL) processes representation, the Noah land surface model (Noah LSM; Chen & Dudhia, 2001) for surface options, and the Kain-Fritsch cumulus parameterization, which was activated in both domains.
Cloud microphysics is treated using the scheme of Morrison et al. (2005), which was further detailed by Morrison et al. (2009) (hereafter M09).The robustness of this scheme has already been evaluated in various atmospheric conditions worldwide, including remote Antarctic (Sotiropoulou, Vignon, et al., 2021) and Arctic clouds (Schäfer et al., 2023), as well as orographic (Billault-Roux et al., 2023;Georgakaki et al., 2022Georgakaki et al., , 2023) ) and more convective cases (Karalis et al., 2022).M09 takes into account two liquid species, namely cloud droplets and raindrops, as well as three ice species: cloud ice, snow, and graupel.All liquid and ice hydrometeors are represented using a double-moment approach, with the exception of cloud droplets, for which a constant number concentration must be prescribed.For this study, we assumed that an average cloud droplet number concentration of 100 cm 3 is suitable for Arctic clouds (e.g., McCoy et al., 2020;Young et al., 2016).Note that our prior work in orographic clouds has shown minor sensitivity of the most important microphysical processes (e.g., riming and SIP rates) and properties (i.e., ice and liquid water content) to the prescribed cloud droplet number (not shown).Regarding ice formation processes, M09 includes parameterizations for homogeneous freezing at temperatures below 40°C, as well as heterogeneous ice nucleation initiated below 4°C.The latter accounts for immersion freezing of cloud droplets and raindrops (Bigg, 1953), contact freezing (Meyers et al., 1992) and condensation/deposition freezing nucleation (Cooper, 1986), all of which are dependent on temperature.
Regarding SIP, there is only one mechanism implemented in the standard version of the M09 scheme-the HM process (Appendix A)-which is activated after cloud droplets or raindrops rime onto snow or graupel particles within a specific temperature range between 8 and 3°C (Reisner et al., 1998).The activation of this process depends on preset thresholds concerning the minimum mass mixing ratios of the involved ice and liquid species that are found to limit its efficiency-especially in polar clouds-and were therefore removed (Atlas et al., 2020;Karalis et al., 2022;Young et al., 2019).
In this study, the M09 scheme was updated to include the description of two additional SIP processes: BR and DS.A detailed description of the BR implementation is provided in Sotiropoulou, Vignon, et al. (2021), while a brief overview of the employed parameterization is also provided in Appendix A. To implement this process, it is important that all possible collisional interactions among the three ice hydrometeor species are described within M09 (see Appendix B of Sotiropoulou, Vignon, et al., 2021).The number of fragments produced is then parameterized as a function of the ice particle size, habit, rimed fraction as well as the collisional kinetic energy, following Phillips et al. (2017).Two types of collisions were considered: those involving only high-density precipitating ice particles (such as graupel), and those involving collisions of rimed snow or cloud ice particles with any ice species.The rimed fraction and ice habit were exclusively used in the formulation describing the latter type of collisions.While M09 does not explicitly predict these two input variables, specific assumptions were made to address this limitation.Starting from the rimed fraction, a value of 0.4 was prescribed, which was found to lead to improved model performance in simulations of polar clouds (Sotiropoulou et al., 2020;Sotiropoulou, Vignon, et al., 2021).To align with this assumed rimed fraction, BR was activated only when there was a nonzero riming rate of either raindrops or cloud droplets onto the ice particle undergoing fragmentation.Concerning ice habit, we assumed planar ice particles due to their ability to capture a broader range of conditions in terms of temperature and particle shapes compared to dendritic ice crystals.Sotiropoulou, Ickes, et al. (2021) suggested that the absence of a thorough handling of ice habit and rimed fraction in most bulk microphysics schemes does not seem to significantly impact the depiction of BR in polar conditions.All secondary ice fragments produced by BR are added to the cloud ice category.
Implementation of the DS mechanism in the M09 scheme of WRF was described in Georgakaki et al. (2022) (summarized here in Appendix A).We considered two collision modes that can cause the freezing and subsequent shattering of large raindrops, as described by Phillips et al. (2018).The first mode, or 'mode 1', occurs when a supercooled raindrop collides with a less massive cloud ice particle or when an INP initiates freezing in immersion mode.Freezing probability is set to unity and zero for temperatures below 6°C and above 3°C, respectively, with intermediate values in between.The shattering probability depends on raindrop size, being 0 for sizes smaller than 50 μm, 1 for sizes larger than 60 μm, and varying for sizes in between.The second mode, or "mode 2," involves collisions between raindrops and larger ice particles such as snow or graupel (James et al., 2021), producing only tiny ice fragments that are passed to the cloud ice category.Big fragments are treated as either graupel, snow, or frozen drop depending on the type of collision that initiated the raindrop freezing process.
The regional climate set-up (Figure 1) with augmented cloud microphysics was employed to run a reference WRF simulation for 2 years from 2016 to 2017, preceded by a 1-week spin-up period starting on 25 December 2015.This time frame was considered sufficient to capture the possible parameter space in the atmospheric state.Threedimensional snapshots of the simulation with a 10-km horizontal grid spacing (the area bounded by the dashed lines in Figure 1) were saved every 3 hr.These snapshots formed the basis for developing the ML parameterization which is described in the following section.

Random Forest (RaFSIP) Parameterization Development
RaFSIP is based on a random forest regressor (RFR) algorithm, which is a supervised ML technique, consisting of an ensemble of decision trees (Breiman, 2001) trained directly from the outputs of the 10-km horizontal grid size simulation.To minimize variance, each decision tree is trained on a random subset of the training data and a random subset of the inputs-also called as features-of that data.The RFR prediction is calculated by averaging over the individual decision trees, which reduces the risk of overfitting and improves generalization performance compared to a single decision tree.RFR can capture non-linear relationships between input features and output variables, but cannot extrapolate due to its predictions being averages over subsamples of the training data set.This helps ensure that an RF-based parameterization will be stable when coupled with a model in an online setting, such as in our case.

Description of the Two Parameterization Approaches
Depending on the temperature and the presence of rime onto ice particles, we expect different SIP processes to be activated.This is the reason why the two versions of the RaFSIP parameterization consist of a number of RFR models, each one of them being responsible for predicting the effect of SIP as dictated by the thermodynamic and microphysical state describing every model timestep.The RFR models comprising the RaFSIPv1 parameterization are trained to predict the IEF, which is the ratio between the production rate of secondary ice particles from each SIP process (HM rate , BR rate and DS rate ) and primary ice particles via heterogeneous freezing of cloud droplets (PIP rate ).Hence, RaFSIPv1 provides 3 predictions: IEF HM , IEF BR , and IEF DS , the contribution of which can be added to yield the total IEF = SIP rate /PIP rate .The notion of IEF is frequently employed in the literature as an indication of the prevalence of SIP, in both observational (e.g., Korolev et al., 2022;Wieder et al., 2022) and modeling studies (e.g., Waman et al., 2022;Zhao et al., 2023).
In RaFSIPv2, the RFR models predict directly the three SIP rates-HM rate , BR rate and DS rate (in kg 1 s 1 ), derived from the 10-km horizontal grid size simulation outputs.The three rates are combined to determine the total SIP rate .As the HM and DS processes move mass from the liquid to the ice species when active (albeit from liquid drops collected by ice in the former case, and frozen liquid drops in the latter), both versions of the RaFSIP parameterization are trained to predict Qc tr and Qr tr (in kg kg 1 s 1 ) corresponding to the mass transferred from cloud droplets and raindrops to the cloud ice category, respectively.Therefore, the two versions of the RaFSIP parameterization can provide up to five different predictions (Figure 1), given certain SIP conditions are met.
Figure 2 provides an overview of the flowchart, illustrating the sequence of steps each RFR model makes for predicting the impact of SIP.We consider two main temperature regimes: 8 ≤ T< 3°C (HM zone) and 25 ≤ T< 8°C, which is essential as the HM mechanism can only be effective within the former one.By combining these temperature ranges with the presence or absence of mass production rate of raindrops collected by ice particles (RIMR in kg kg 1 s 1 ; Equation 1), we obtain 4 different RFR models employed by both versions of the RaFSIP parameterization (Tables 1 and 2).Two models can be activated inside the HM zone, namely "forestALL" and "forestBRHM", that are replaced by the "forestBRDS" and "forestBR" at lower temperatures (the suffix of the model name denotes the SIP processes in action).The final selection of which forest to use in the corresponding temperature range is based on the presence of a nonzero RIMR at the current timestep, which is the criterion used to identify the cases where DS might be taking place.The two forests that take into account the effect of DS are the forestALL and forestBRDS models.
Both RaFSIP parameterizations incorporate SIP through BR and DS at temperatures as low as 25°C, which is supported by recent observational findings (Korolev et al., 2022;Pasquier et al., 2022;Wieder et al., 2022).Approaching higher subzero temperatures (T > 4°C), RaFSIPv1 does not predict ice enhancement, since heterogeneous freezing of cloud droplets (i.e., the denominator in the IEF expression) does not occur in the M09 scheme (Section 2.1).In contrast, RaFSIPv2 can still account for the effect of SIP as it is not directly linked to the PIP rate .An additional RFR model called "forestBRwarm" is therefore used in RaFSIPv2 (inset of Figure 2), to calculate the BR rate at temperatures between 3 ≤ T ≤ 0°C.At these warm subzero temperatures, DS is also proposed to contribute to SIP through recirculation (Korolev et al., 2020), but only in tropical or midlatitude frontal systems (Lauber et al., 2021) owing to strong convective turbulent updrafts that are rare in polar stratiform clouds.Therefore, RaFSIPv2 considers that only the BR mechanism initiated by collisions of seeding rimed ice particles can be efficient at temperatures above 3°C.

Features and Training Data Set
Both RaFSIP versions share the same input features, including ambient temperature (T in K), total ice water content (IWC in kg kg 1 ), liquid water content (LWC in kg kg 1 ), relative humidity with respect to ice (RHI) and mass production rate of cloud droplets rimed onto ice particles (RIMC in kg kg 1 s 1 ; Equation 2), which are among the features that have been used in all RFR models of the two versions of the parameterization.Since the freezing of raindrops leading to SIP through DS is restricted to instances with a nonzero mass of raindrops collected by ice species, RIMR (Section 2.3) is included as an input in models accounting for its effect (i.e., forestALL and forestBRDS).With these six inputs we seek to effectively capture all range of cloud states that can initiate and drive ice multiplication.Despite the prescribed cloud droplet number in M09, we anticipate the two riming rates, RIMC and RIMR, to serve as more relevant indicators of the potential for SIP.Notably, the cloud droplet number is not considered an input in the inline prognostic SIP parameterizations for predicting SIP fragments, in contrast to the riming rates.The training data set of RaFSIP was constrained to cases with substantial ice and liquid water content varying between 10 5 and ∼1 gm 3 , accompanied by nonzero riming tendencies, to ensure the effectiveness of SIP.The range of both riming rates within the 2-year training data set was between 10 18 and 10 6 kg kg 1 s 1 .The 6 inputs are chosen in a way to simplify the detailed SIP formulations (e.g., Phillips et al., 2017Phillips et al., , 2018), yet are expected to generate complex mappings for training the RFR models.Here we use the instantaneous predictions of the model for all individual vertical grid cells as inputs for all RFRs, rather than their vertical profiles (e.g., O'Gorman & Dwyer, 2018;Yuval & O'Gorman, 2020).We generated the training data set for all RFR models using the 3-hourly WRF outputs with a horizontal grid spacing of 10-km (Section 2.1), extracted from four 2 × 2 grid cell regions (16 grid cells in total) indicated by purple triangles in Figure 1.Half of these grid cells represented continental conditions, while the other half were over the sea.The selection of these 4 particular locations was arbitrary, but the maximum number of training samples was defined after hyperparameter tuning (Section 2.2.3).The absence of a latitudinal preference in SIP rates between the Northern and Southern hemispheres (Zhao & Liu, 2021) suggests that a limited number of grid points from the 2-year finer-resolution simulation of WRF is adequate for capturing the thermodynamic and microphysical conditions associated with SIP events in stratiform conditions.
Horizontal and vertical coarse-graining is frequently applied when developing ML parameterizations for GCMs based on high-resolution simulations (e.g., Brenowitz & Bretherton, 2019;Grundner et al., 2022;Yuval & O'Gorman, 2020).This scale-chain approach is essential when ML algorithms are trained on high-resolution models with convection-permitting resolutions but are intended for implementation in coarser-scale models where convection is parameterized (e.g., Yuval & O'Gorman, 2020).A key difference between RaFSIP and ML algorithms designed to speed up convection parameterizations in GCMs is that it has been trained on a data set derived from a 10-km grid spacing domain, where convection is already parameterized using a cumulus scheme.Therefore, our goal is to omit coarse-graining and investigate whether a SIP parameterization derived from a 10km horizontal grid size domain (innermost WRF domain) remains effective at the coarser grid spacing of 50-km in the parent domain of WRF, which aligns more closely with the resolutions seen in GCMs.We acknowledge that the scale dependency of some inputs may introduce uncertainties to RaFSIP predictions.However, we observed consistent riming rates and mass concentrations, whether derived from the 10-km grid spacing domain or the coarse-grained 50-km one (not shown).Among the input variables, LWC exhibited the highest resolution dependency, with values that can be up to one order of magnitude lower when derived from the 50-km grid size domain compared to the 10-km one.This, however, impacted only predictions related to HM, with the remaining SIP processes remaining insensitive to such scale differences.Considering that HM operates within a specific temperature range, and its contribution may be overshadowed by a more efficient BR mechanism (Section 3), we anticipate that the overall performance of RaFSIP, especially under stratiform conditions, will not be significantly impacted when data from a coarser grid spacing domain than the one RaFSIP has been trained on are used as inputs.Further evaluation of the performance of RaFSIP with coarser resolution inputs is described in Section 4.
To prepare for training, each feature and output variable is normalized using the natural logarithm, which is important given the range of several orders of magnitude they span (Figure 3).During the training process, we excluded IEF values below 10 1 and SIP rates below 10 5 kg 1 s 1 (Figure 3).This was made because such values would not result in significant SIP, and we aimed to reduce the range of predictions.The training data set comprises 20 months of 10-km horizontal grid spacing WRF outputs, representing approximately 83% of the data.
The remaining 17% of the data, consisting of 4 months (January, April, July, and October 2017), were used for testing the offline performance of the RaFSIP parameterization (Section 3).The RandomForestRegressor class from the scikit-learn package version 1.2.0 (Pedregosa et al., 2011) was used to train each RFR model.All RFR models were stored as ASCII files.

Choice of Hyperparameters
To improve the performance of the RFR models, different hyperparameters governing the learning process can be tuned.The Mean Squared Error (MSE) resulting from the RFR predictions was evaluated using 10-fold cross validation on the training data set.We selected the number of trees in each forest, the minimum samples per leaf node, and the number of training samples as the three most important hyperparameters to tune. Figure S1 in Supporting Information S1 shows the validation curves for the forestALL model (used in RaFSIPv2), which illustrate how the performance of the model varies with different hyperparameter values.The final decision on hyperparameters involves a trade-off between model complexity and runtime, with the latter being crucial for optimizing online performance when coupled with the WRF model.We chose 10 trees in each RFR, a minimum of 4 samples per leaf node, and 25,000 training samples for all RFR models, as these produced comparable validation curves.Note that, for online applications of the RFR models, it is common practice to use 10 decision trees to ensure computational efficiency (O'Gorman & Dwyer, 2018;Yuval & O'Gorman, 2020).Focusing on the number of training samples, Figure S1 reveals that the MSE of our predictions plateaued beyond approximately 25,000 training samples.Instead of randomly shuffling and selecting 25,000 samples, we deliberately chose the 16 grid cells shown in Figure 1, ensuring that the ML algorithm was trained across various locations, covering different seasons throughout the 2-year simulation period.The offline performance of the RaFSIP parameterizations will be discussed in Section 3, followed by the online performance in Section 4.

Implementation of RaFSIP Parameterization in WRF
To implement the RaFSIP parameterization in WRF, a Fortran 90 module is developed to read and store the parameters for building the RFR models.The ASCII files are only read during the first model timestep and all forest parameters are passed as public variables into the microphysics code.Within the M09 microphysics routine, the RaFSIP parameterization replaces the three detailed inline SIP parameterizations (Section 2.1), taking the fields of T, RHI, LWC, IWC as well as RIMC and RIMR as inputs to make predictions.Calculations for LWC and IWC, consider contributions from all liquid (cloud droplet and raindrops) and ice (cloud ice, snow and graupel) species.RIMC and RIMR account for the various ways cloud droplets or raindrops are collected by frozen particles.Specifically, within the M09 microphysics scheme, RIMR is defined as follows: where each term represents the bulk rate (in kg kg 1 s 1 ) at which the raindrop mass is collected by snow (rs), graupel (rg) or cloud ice (ri) particles.Similarly, RIMC is defined as: where again the three terms correspond to the three mass tendencies (in kg kg 1 s 1 ) of cloud droplets accreted on snow (cs), graupel (cg) or cloud ice (ci).At each model timestep and grid cell, if certain requirements are met (discussed in Section 2.2.1), a different RFR model may be activated to represent the effect of single or combined SIP processes.
At the end of each model timestep, the RaFSIP parameterization calculates the total SIP rate (in kg 1 s 1 ) either directly or through the use of the IEF, depending on the chosen version.In RaFSIPv1, SIP rate is obtained by multiplying the number tendency of heterogeneously frozen cloud droplets (in kg 1 s 1 ) by the total IEF.However, this approach has a significant caveat: SIP cannot be accounted for unless a nonzero PIP rate is predicted by the model.Therefore, RaFSIPv1 can only be used at temperatures colder than 4°C, where PIP is initiated in the model.To ensure that the final SIP rate does not exceed a certain threshold, an upper limit of 100 particles kg 1 s 1 was applied, corresponding to the highest simulated value in the training data set.While our results were not sensitive to the choice of this upper limit, it is important to note that, as RaFSIPv1 is connected with the PIP rates produced by M09, employing this upper limit is recommended when applying RaFSIPv1 in other microphysical schemes.The SIP rate predicted by RaFSIPv1 is then added to the cloud ice number conservation equation.The RaFSIPv2 parameterization follows a similar approach to RaFSIPv1, but does not link SIP rate to the instantaneous PIP rate , rendering the scheme less dependent on the specific microphysics scheme on which it was trained.This allows RaFSIPv2 to be applicable over a wider temperature range from 0°C to 25°C, without requiring a nonzero PIP rate .While RaFSIPv1 exhibits more restricted applicability in the model compared to RaFSIPv2, the use of IEF remains promising for quantifying and parameterizing the effect of SIP in MPCs, as it holds potential for direct comparisons with field observations, a capability not provided by the SIP rates predicted in RaFSIPv2.
When HM and/or DS mechanisms are active, the transported masses Qc tr and Qr tr are subtracted from their respective cloud categories and added to the conservation equation for cloud ice mass mixing ratio.It is important to note that, in the absence of SIP, these masses would have been transported to the mass conservation equation of the corresponding rimed ice particles.For instance, the first term in Equation 2 represents the mass transfer from cloud droplets to snow after riming occurs.However, when SIP is active, a portion of this rimed mass (Qc tr ) is redirected to the small cloud ice fragments.To avoid double-counting of the transferred liquid masses, Qc tr and Qr tr must be deducted from the riming tendencies of the involved ice species.To address this, in cases where more than one term of Equations 1 and 2 is positive, we subtract the transferred masses from the term with the greatest contribution.Considering that the terms involved in RIMC and RIMR (Equations 1 and 2) can vary by several orders of magnitude, this approach ensures that the simulations will not be subject to significant biases.
To evaluate the online performance of the new RaFSIP parameterizations when coupled with the WRF model, we conducted a 1-year simulation using the same set-up depicted in Figure 1 (Section 2.1).This simulation covers the period between September 2019 and August 2020, which is more than 1.5 years after the training period, with one additional week of spin-up starting from 25 August 2019.The testing simulation was conducted three times.First, we carried out the "CONTROL" simulation, which does not take into account any SIP process.Then, we performed the "ALLSIP" simulation, which includes all the detailed inline parameterizations of SIP outlined in Section 2.1.Finally, the "RaFSIP" simulation refers to the one where the RaFSIP parameterizations replaced the detailed SIP descriptions.By comparing the ALLSIP and RaFSIP simulations in Section 4, we seek to examine the robustness of the new SIP parameterizations in terms of computational efficiency and simulation quality.

Offline Performance of RaFSIP Parameterizations
To evaluate the performance of both RaFSIP parameterizations, each RFR model was tested on the 4-month data set that was not used during training (Section 2.2.2).The performance was assessed using two metrics summarized in Tables 1 and 2: Root-Mean-Square Error (RMSE) and coefficient of determination (R 2 ), which were calculated based on the normalized output variables.R 2 is defined as one minus the ratio of the sum of squared errors to the true variance.The offline performance of the forestALL model, which is used by both RaFSIP parameterizations when all conditions for SIP processes are met, is shown in Figure 3.For the rest of the RFR models, the offline performance is displayed in Figures S2 to S5 in Supporting Information S1.
All five outputs of the forestALL model present a wide range of values (Figure 3), making accurate prediction challenging.For instance, IEF BR and IEF HM vary from 10 1 to 10 11 , while IEF DS is confined below ∼10 8 , with predictions more spread out along the 1:1 line, especially at lower values.The RaFSIPv1 parameterization predicts an inverse relationship between IEF and PIP rate (not shown), primarily due to the PIP rate being in the denominator of the IEF expression.The competing relationship between SIP and PIP has also been highlighted in previous modeling studies (e.g., Georgakaki et al., 2022;Zhao & Liu, 2022).Note that the predicted IEF values are higher than those frequently reported in the literature (e.g., Waman et al., 2022;Zhao et al., 2023), as they are extracted from the ratio of tendencies rather than integrated particle concentrations.Turning to the offline predictions of SIP rates from RaFSIPv2 (Figures 3f-3h), we observe a distribution around the one-to-one line that aligns with the IEF predictions from RaFSIPv1.Despite methodological disparities between the two approaches, RaFSIPv2 predictions underscore once again the greater efficiency of BR and HM in generating SIP particles compared to DS. Notably, the maximum production rates of DS remain at least one order of magnitude lower than those achieved by the other two processes.
Although RaFSIPv1 and RaFSIPv2 use the same predictors, the former is limited to cases where PIP rate is nonzero, resulting in noticeable differences among the two approaches.To evaluate the relevance of each input feature in the predictions made by the RFR models, we employed the permutation importance metric (Breiman, 2001;McGovern et al., 2019) from the RandomForestRegressor class of scikit-learn (Figures S6 to S10 in Supporting Information S1).This metric estimates feature importance by randomly permuting its values and measuring the consequent decrease in model performance.This approach can capture the effect of both linear and nonlinear relationships between features and target variables, and can also help identify any potential overfitting issues in the RFR models.Note that here we used the negative RMSE as the scoring parameter and normalized the resulting variable importance by dividing by the maximum value among features.
Examining the permutation feature importance reveals that the most significant features in the predictions across all RFR models are the riming rates and IWC.Here it is important to highlight a potential limitation of the employed feature importance measure, as resulting scores may be less reliable when handling highly correlated input features.This arises from the fact that permuting one correlated feature might not substantially affect the performance of the model if the other compensates for its importance.In our case, the correlation matrix of the inputs (not shown) revealed a strong positive correlation (correlation coefficient of 0.7) between the two riming rates, RIMC and RIMR.This correlation may influence the interpretability of permutation importance, especially for the forestALL (Figure S6 in Supporting Information S1) and forestBRDS (Figure S8 in Supporting Information S1) models, which utilize both variables as input features.Inside the HM zone, RIMC emerges with the highest score for the forestALL (Figure S6 in Supporting Information S1) and forestBRHM (Figure S7 in Supporting Information S1) models employed in both RaFSIP versions, which is not surprising given that the two riming rates are explicitly used as inputs in the detailed parameterized expression of HM in the M09 scheme (Section 2.1).Outside the HM zone, where BR dominates the SIP rates compared to DS, IWC becomes the most important feature in the predictions of the RaFSIPv2 models, which is particularly evident for forestBR (Figure S9 in Supporting Information S1) and forestBRwarm (Figure S10 in Supporting Information S1).An increase in IWC, associated with ice sedimenting from higher-level clouds, is frequently identified as a driving force for the BR process, as supported by existing literature (Georgakaki et al., 2022;Järvinen et al., 2022;Pasquier et al., 2022;Ramelli et al., 2021;Sotiropoulou, Vignon, et al., 2021).
Using the different RFR models depending on the prevailing thermodynamic and microphysical conditions allowed us to evaluate the performance of the RaFSIP parameterizations under various scenarios.The distribution of all RFR predictions against true values, as shown in Figure 3 and Figures S2 to S5 in Supporting Information S1, indicates that predictions with a higher relative frequency of occurrence are closer to the 1:1 line for both RaFSIP approaches.The statistical metrics presented in Tables 1 and 2 can differ significantly depending on the prevailing conditions and the SIP processes involved.The RMSE for RaFSIPv1 predictions varied from ∼10%-45%, while for RaFSIPv2, it ranged between ∼15%-40%.The transferred masses due to SIP, Qc tr and Qr tr , are generally well-predicted in most of the RFR models with an RMSE of less than 30%, despite spanning more than 10 orders of magnitude.
In cases where multiple SIP processes are present, the prediction error tends to be higher for DS, as shown in Figures 2c and 2h and reflected in the higher RMSE values in Tables 1 and 2. To improve the offline accuracy of the forestALL and forestBRDS models for predicting this SIP mechanism, we tested the inclusion of the size of raindrops undergoing DS as a potentially important predictor.The resulting RMSE metrics showed a ∼5% decrease in the forestBRDS model, while the forestALL model remained unaffected (not shown).However, since the total effect of SIP is cumulative and DS is the least significant mechanism compared to the other two, any errors associated with its uncertain representation are not likely to significantly impact the online performance of the RaFSIP parameterizations.Therefore, we argue that adding more predictors may not be worth the increased model complexity.
When only BR is active, either in the colder (i.e., forestBR) or warmer (i.e., forestBRwarm) regime, the RMSE scores range between ∼30%-45% (Tables 1 and 2).This could indicate that an important input feature related to BR may be missing or that the relationship between the input features and the target variable (IEF BR or BR rate ) is inherently difficult to model accurately due to the level of non-linearity.In the physically based parameterization by Phillips et al. (2017), the SIP fragments generated after BR depend on collisional kinetic energy, which is an uncertain parameter as it is a function of the difference in terminal velocities of colliding ice particles.The latter follows a fall-speed-diameter relationship, the parameters of which are highly uncertain and can directly impact both the precipitation rates and the SIP efficiency (Karalis et al., 2022).Therefore, we chose not to include collisional kinetic energy in our set of predictors, as it would increase model complexity without necessarily improving accuracy.The predictors we chose strike a good balance between model complexity and accuracy, allowing us to describe SIP occurrence as a function of simplified inputs that most GCMs can predict.

Coupling the RaFSIP Parameterization With WRF
In this section, we aim to assess the performance of the parameterization when coupled back with the dynamical core of the WRF model, making fast and reasonable predictions at runtime.To facilitate the discussion, in the rest of the paper we will focus on the RaFSIPv2 simulation, while the corresponding RaFSIPv1 simulation outputs will be available in the Supporting Information.RaFSIPv2 can predict SIP rates directly, making it less coupled with the WRF model on which it was trained.In contrast, the approach followed in RaFSIPv1 tends to be more dependent on the PIP rates predicted by each model, implying that its good online performance in WRF may not be guaranteed in other models using different PIP schemes.
When introducing RaFSIP into the M09 scheme of WRF, two main concerns had to be addressed in order to ensure the robustness of the approach: model stability and computational efficiency.The stability of the model was tested during the 1-year test simulation starting from September 2019 until August 2020 (Section 2.1).The coupled RaFSIP-WRF model operated smoothly and consistently across diverse conditions, without experiencing model blow-ups.Regarding computational efficiency, the new parameterization performed comparably to the inline SIP parameterizations without increasing computational demands.During the same simulation period, both ALLSIP and RaFSIP required approximately 10% more elapsed real-time compared to CONTROL to complete.This implies that the RaFSIP scheme can be easily incorporated into existing models without significant impact on runtime.In the following discussion all results will be divided into four groups according to season.

Horizontal Distribution of Ice Crystal Number Concentrations
Latitude-longitude snapshots of the median ICNCs obtained from the instantaneous 3-hourly outputs of the CONTROL, ALLSIP, and RaFSIP experiments (Section 2.3), along with the resulting R 2 between the predictions of the latter two, provide a way to assess the quality of WRF predictions when coupled with the RaFSIPv2 parameterization (Figure 4).Note that the modeled ICNCs include the contribution of all three ice species (cloud ice, snow and graupel).The horizontal distribution of ICNCs depicted in Figure 4 was determined by calculating median values across both the time and altitude dimensions.For this calculation we considered only the in-cloud conditions represented by ICNCs >10 5 L 1 and temperatures between 25 and 0°C, which is the range where RaFSIPv2 can be active.The subsequent calculation of R 2 scores was performed for each latitude and longitude, after averaging the RaFSIPv2 and ALLSIP simulation outputs over time.To facilitate the identification of horizontal patches where RaFSIPv2 demonstrates enhanced or decreased predictive skills in terms of R 2 , we further aggregated the information from the coarser-resolution grid mesh by dividing the 50-km WRF domain into non-overlapping 3 × 3 sub-grids, and computing the average value of each sub-grid.Thus the R 2 map in Figure 4 is projected onto a 49 × 49 mesh.The RaFSIPv2 parameterization coupled with WRF can skillfully predict the median horizontal distribution of ICNCs, closely following the pattern produced by the ALLSIP simulation with the detailed SIP descriptions.Darker blue shades are produced by ALLSIP and RaFSIP indicating the higher median ICNCs compared to CONTROL.The increase in median ICNCs is observed for all four seasons, even though it is more prevalent during winter and spring, particularly in the sea region between Greenland and Canada to the west and extending to the north of the Scandinavian Peninsula to the east.In these areas the median seasonal ICNCs can exceed 10 L 1 , suggesting the significant role of ice multiplication processes in elevating the median ICNCs by up to an order of magnitude when SIP conditions are met.
Despite being trained on a small subset of 10-km horizontal grid spacing WRF outputs, the RaFSIPv2 parameterization can generalize and accurately predict the locations of ice enhancement in the coarser grid spacing domain as dictated by the ALLSIP simulation.The R 2 maps superimposed in Figure 4 reveal patches of especially high predictive skill denoted with R 2 scores exceeding 90%.The accuracy of the predictions largely exceeds 80% of R 2 over localized areas of significant ice enhancement, shown with darker shades on Figures 4g and 4k.This implies the ability of RaFSIPv2 to reproduce the physical processes occurring in polar MPCs.
Some localized areas with lower R 2 scores are primarily observed near the edges of the coarser WRF domain, particularly at latitudes below 60°N (e.g., Figure 4l).One important factor affecting these patterns might be the geographical distribution of the RaFSIP training data set, with the majority of grid points selected for training extracted from latitudes higher than 70°N.Consequently, the lower performance in the midlatitudes (<60°N) may be anticipated due to the Arctic-centric focus of the training data set.The southeast part of Greenland also presents challenges for accurate model predictions, leading to R 2 dropping below 60% during winter, in a region where RaFSIP tends to underestimate ICNCs compared to ALLSIP (not shown).In this region, enhanced winter precipitation is associated with frequent cyclogenesis and the presence of the Icelandic low.The suboptimal performance of RaFSIPv2 in this specific area may be attributed to the underestimation of SIP rates, particularly in connection with more convective systems.
Compared to other seasons, the predictive skill of RaFSIPv2 is reduced during summer, with R 2 dropping below 40% in very localized regions close to the edges of the WRF domain.Higher temperatures are not likely the cause of lower R 2 values, as the RaFSIPv2 models were trained using a wide range of temperatures, and temperature was not found to be among the most important input features for the predictions of the RFR models (Figures S6-S10 in Supporting Information S1).A possible explanation could be that during summer, the expected shift toward warmer clouds may lead to increased riming rates, which can be associated with higher SIP rates in RaFSIPv2.Replacing the RaFSIPv2 parameterization with RaFSIPv1 in WRF results in a further decline in predictive skill in all four seasons, as shown in Figure S11 in Supporting Information S1.Note that the R 2 metrics are now derived with a focus on temperatures between 25 and 4°C, which is the range where PIP rate can be nonzero allowing RaFSIPv1 to be active.The R 2 scores decrease to below 60% in most regions with significant ice enhancement, with a further drop to below 40% during the summer months.The degradation in the performance of the model is most likely because RaFSIPv1 requires a positive PIP rate to provide ice enhancement.
Figure 5 presents median statistics extracted from all three sensitivity experiments during the cold (autumn and winter) and warm (spring and summer) seasons.The latitudinal zones north of 70°N are divided into three groups, and the results are grouped into two distinct temperature ranges ( 10 ≤ T ≤ 0°C and 20 ≤ T < 10°C) where SIP effects are more prevalent.During the cold period, it is observed that the CONTROL simulation leads to higher median ICNCs in the warmer temperature range compared to the colder one at latitudes 70-80°N (black circles in Figures 5a and 5c), which may seem contradictory since PIP efficiency generally increases with decreasing temperatures.This could potentially be explained by ice seeding from colder clouds enhancing the ICNCs in lower-level warmer clouds at temperatures above 10°C, or by increased efficiency of ice aggregation leading to significantly lower median ICNCs in the colder range ( 20 ≤ T < 10°C) (Barrett et al., 2019;Chellini et al., 2022).
Introduction of SIP in WRF results in a substantial increase in median ICNCs, particularly in the higher temperature bin, with ALLSIP and RaFSIP simulations showing ICNCs ∼2-3 and ∼3-6 times higher than the CONTROL simulation for cold (Figures 5a and 5c) and warm (Figures 5b and 5d) periods, respectively.On the other hand, in the lower temperature bin, the ICNC enhancement in ALLSIP and RaFSIP is relatively modest, reaching only up to ∼20%-25% during both examined periods.It is worth noting that latitudes above 80°N display systematically lower ICNCs, with no significant contribution of SIP being observed at temperatures below 10°C, likely due to complete cloud glaciation.Overall, median statistics derived from the WRF simulation coupled with RaFSIPv2 are in good agreement with the ALLSIP simulation during the cold and warm periods examined here (Figure 5).
To assess the agreement between our WRF simulations and satellite remote sensing data, we compare our median values from the 3 sensitivity experiments to the 10-year median values between 2006 and 2016 presented in Papakonstantinou-Presvelou et al. (2022) (hereafter denoted as PP22).Their study focused on low-level pure ice clouds and demonstrated that higher ICNCs are found over sea ice compared to open ocean-a finding contrary to previous expectations.Their observations are superimposed on Figure 5 with gray markers.To enhance comparability with these observations, we calculated the median ICNC statistics focusing on grid points over sea ice with a glaciation fraction exceeding 98% and a sea ice concentration higher than 50%.These are indicated by the filled rhombuses on Figure 5.The glaciation fraction was determined by calculating the ratio of modeled IWC to total water content (TWC = IWC + LWC).While our simulations fall outside the 10-year medians extracted from PP22, these observations are still valuable for a qualitative analysis.
In the glaciated clouds simulated between 10°C and 0°C, the inclusion of SIP in the ALLSIP and RaFSIP schemes brings the modeled ICNCs closer to the 10-year median observed values.At these higher subzero temperatures, both ALLSIP and RaFSIP enhance the median ICNCs by a factor of up to ∼2 compared to CONTROL during both periods (Figure 5).In glaciated conditions, the efficiency of SIP is constrained, leading to less pronounced ice enhancement.This is particularly evident at temperatures below where the efficiency of ALLSIP and RaFSIP in producing SIP particles diminishes.Under these temperature conditions, all three sensitivity experiments tend to underestimate the observed median ICNC levels, which can reach up to ∼3 or 4 L 1 at all latitudes.In moderately cold clouds, where SIP is expected to dominate total ice production rates over PIP, especially in single-or multi-layer Arctic clouds (Zhao & Liu, 2022), uncertainties in SIP representation may contribute significantly to the simulated ICNC biases.Additionally, biases associated with the representation of PIP in the model cannot be overlooked.

Joint Probability Distributions of Simulated Cloud Properties
Figure 6 illustrates the bivariate joint probability density function (PDF) of the median ICNCs as a function of temperature and glaciation fraction, which is a normalized quantity that can effectively reveal how much of the cloud mass is ice.This figure aims to provide further insight into the conditions where the RaFSIP parameterization is expected to perform well compared to the ALLSIP simulation.The statistics were derived using instantaneous 3-hourly outputs obtained from the 50-km grid spacing domain of the CONTROL, ALLSIP, and RaFSIP experiments (Section 2.3), focusing again on the in-cloud conditions characterized by significant concentrations of ice crystals (ICNC >10 5 L 1 ), as well as ice and liquid masses (IWC, LWC >10 6 gm 3 ).

Median Ice Crystal Number Concentration
Figure 6 reveals that the WRF model, coupled with the RaFSIPv2 parameterization, can accurately capture the most prominent patterns of the median 2D-binned ICNCs in all four seasons compared to the ALLSIP simulation with detailed inline microphysics.This is particularly true for temperatures ranging from 8°C to 3°C and glaciation fractions below ∼25%-30% or higher than ∼80%.The turquoise bands in the lower right of the 2Dbinned ICNCs indicate the importance of SIP in relatively warm polar clouds.At these temperatures, supercooled liquid water prevalent in simulated MPCs enhances riming tendencies (RIMC and/or RIMR) and facilitates the action of SIP, resulting in modeled ICNCs up to an order of magnitude higher in the ALLSIP and RaFSIP simulations compared to CONTROL.These findings align with the simulations of Arctic clouds presented by Sotiropoulou et al. (2020), which reported a 10-20 fold increase in ice when both HM and BR were considered, with the effectiveness of BR decreasing in moderately colder clouds (Sotiropoulou, Ickes, et al., 2021).
At higher glaciation fractions (between ∼25%-80%) within the same temperature range, we observe that the RaFSIPv2 simulation tends to overestimate the median ICNCs predicted by ALLSIP.Indeed, from the normalized histograms of 2D-binned median ICNCs superimposed in Figure 6, we can also infer that RaFSIPv2 predictions can be up to a factor of 3 higher than the ALLSIP predictions, especially in the lower ICNC ranges (e.g., Figures 6d and 6h).In these temperature conditions, only the forestALL and forestBRHM can be activated, with the latter contributing presumably more due to the expected limited presence of large raindrops, which prevents the frequent activation of the forestALL model.The permutation importance of forestBRHM (Figure S7 in Supporting Information S1) suggests that RIMC and IWC are the two most important features for the predictions of this model.If higher SIP rates are associated with higher IWC and, hence, glaciation fraction, this could explain the slightly overestimated RaFSIPv2 predictions.
At temperatures below ∼ 15°C, where the relative contribution of PIP starts becoming more important than SIP (Zhao et al., 2023), all three sensitivity experiments produce similar ICNC patterns (Figure 6).The vertical zones of constant median ICNCs observed in all three sensitivity simulations originate from the PIP scheme of WRF being dependent only on temperature (Section 2.1).Nevertheless, reduced ICNCs are simulated during the summer months in ALLSIP (Figure 6n) and RaFSIP (Figure 6o) simulations compared to CONTROL, for clouds formed at temperatures below 18°C and glaciation fractions higher than 40% (upper left part of Figures 6n and  6o).At such low temperatures, the competition between SIP and PIP may limit the amount of water available or the distribution of ice crystals throughout the cloud, thus limiting the ice crystals produced by SIP.
Focusing on the upper right corner of Figure 6, where the temperatures are between 5 and 0°C and the glaciation fraction exceeds 80%, we can observe higher ICNCs simulated by the CONTROL simulation than the two experiments that include SIP parameterizations, especially during winter and spring.At temperatures higher than 4°C, elevated ICNCs in the CONTROL experiment can only be explained by ice sedimentation, where precipitating ice particles fall from higher model grid cells.The implementation of SIP in the model can induce numerous ice crystals-albeit with smaller sizes-in the overlying colder clouds or colder parts of the same cloud (Georgakaki et al., 2022;Sotiropoulou et al., 2024).This shift in the ice particle size distribution leads to smaller ice crystals leading to the reduced sedimentation rate of frozen hydrometeors observed in the ALLSIP and RaFSIP simulations.This is further supported by the joint PDF of the modeled IWC values (Figure S12 in Supporting Information S1), where both simulations accounting for SIP result in comparable or even lower IWC values compared to the CONTROL simulation.This can have important implications for the ice and liquid phase partitioning in polar MPCs, determining their microphysical evolution and radiative properties (Curry et al., 1996;Shaw et al., 2022;Tan & Storelvmo, 2019).Overall, the RaFSIP simulation reproduces the most important features in the 2D-binned ICNCs (Figure 6) and IWC (Figure S12 in Supporting Information S1), aligning closely with the predictions of the ALLSIP simulation within a factor of 3.This suggests that the newly developed ML parameterization can qualitatively reproduce the underlying physics of the simulation with the detailed SIP microphysics.
Despite the lower R 2 scores in the horizontally averaged ICNCs (Figure S11 in Supporting Information S1), replacing the RaFSIPv2 parameterization with RaFSIPv1 in the model leads to comparable results when plotted in the joint temperature-glaciation fraction spectrum (Figure S13 in Supporting Information S1).However, the distribution around the 1:1 line is more spread out, especially in spring and summer (Figure S13j and S11p in Supporting Information S1 vs. Figures 6j and 6p).The main difference between the two approaches is at relatively high subzero temperatures (T > 4°C), where RaFSIPv1 disregards completely the effect of SIP due to the zero PIP rates.As a result, the upper right part of the RaFSIPv1 ICNC distribution is more comparable to the one produced by the CONTROL simulation rather than the ALLSIP simulation (Figure S13 in Supporting Information S1).

Median Cloud Liquid Water Content
SIP often causes changes in modeled liquid cloud droplets due to increased number concentrations of ice particles (e.g., Georgakaki et al., 2022;Zhao & Liu, 2021).In addition to direct liquid-to-ice mass transitions caused through the HM and/or DS mechanisms, the elevated concentrations of secondarily formed ice particles can further deplete the surrounding liquid phase hydrometeors in MPCs after their initial growth through vapor deposition, through the action of the Wegener-Bergeron-Findeisen (WBF) or the riming processes.It is therefore essential to investigate how SIP affects the joint PDF of LWC in the ALLSIP simulation and whether the RaFSIP simulation can reproduce the results.
From the complex patterns of LWC produced (Figure 7), we can infer the non-linearity of the microphysical interactions occurring in mixed-phase and ice clouds.The competition between liquid and ice is highlighted by the decreasing cloud LWC as the glaciation fraction increases at a given temperature, which is further anticipated given their interdependence.The dark blue shades confined in the lower part of Figure 7 indicate a persistent supercooled liquid layer with median LWC values exceeding 10 2 gm 3 .This layer is present in all three sensitivity simulations at temperatures below ∼ 4°C with a degree of glaciation below 20% in all seasons except summer, where in the CONTROL simulation it extends up to almost 40% (Figure 7m).This remark can be consistent with the typical structure of high-latitude MPCs, comprising single or multiple stratiform layers of supercooled liquid water from which ice crystals form and precipitate (e.g., Morrison et al., 2012;Shupe, 2011).
Including SIP in the WRF simulations slightly restricts the extent of the liquid layer to lower glaciation fractions in both ALLSIP and RaFSIP, most notably during the summer months.At temperatures above 15°C and glaciation fractions below 60% polar MPCs are produced with reduced LWC compared to the CONTROL simulation.In CONTROL, for instance, during fall and spring, a drop in the simulated LWC to below ∼10 4 gm 3 occurs at temperatures around 8°C only when glaciation fraction exceeds 60% (Figures 7a and 7i).However, in the ALLSIP (Figures 7b and 7j) and RaFSIP (Figures 7c and 7k) simulations, such low LWC values can be attained in clouds with a degree of glaciation as low as 40%.The regions of significant ice enhancement in the ALLSIP and RaFSIP simulations (i.e., turquoise bands in Figure 6) do not overlap with the regions of decreased LWC in Figure 7.This can be explained by the time needed for small secondary ice particles to grow through vapor diffusion before becoming large enough to affect LWC.As glaciation fraction increases, they can gain mass more efficiently at the expense of the surrounding evaporating cloud droplets (WBF) or through their collection and subsequent freezing (i.e., riming).Both simulations accounting for SIP exhibit consistent agreement in the location of the most significant ice enhancement (Figure 6), followed by depletion of the mass of supercooled water (Figure 7).In the region characterized by elevated ICNCs resulting from more frequent ice seeding events (i.e., upper right part of Figure 6a), the CONTROL simulation predicts higher LWC values compared to ALLSIP and RaFSIP.This could be attributed to the presence of numerous ice crystals with higher IWC values (Figures S12a, S12e, S12i, S12m in Supporting Information S1), which sediment rapidly, thereby limiting their ability to compete with supercooled liquid droplets for the available water vapor.In this region, the ALLSIP and RaFSIP simulations predict lower mass mixing ratios for both liquid (Figure 7) and ice hydrometeors (Figure S12 in Supporting Information S1).
At temperatures below approximately 15°C, the modeled patterns of LWC appear consistent across all simulations.At these temperatures the decrease in LWC with increasing glaciation fraction is likely attributed to enhanced PIP rates, such as immersion or contact freezing, even though deposition freezing nucleation is expected to dominate at such low temperatures (not shown).The combination of heterogeneous freezing, WBF, riming and SIP likely controls the complex LWC patterns observed between 20°C and 4°C.Despite the complexity of the simulated system, the RaFSIP predictions align closely with the 1:1 line in the produced normalized histograms, particularly during fall and winter (Figures 7d and 7h), indicating excellent agreement with the ALLSIP simulation.These findings suggest that the RaFSIP parameterization can accurately capture the cloud phase partitioning produced by the ALLSIP simulation, with an uncertainty of up to a factor of 3. The LWC outputs from the RaFSIPv1 simulation (Figure S14 in Supporting Information S1) demonstrate that the use of the IEF for parameterizing the effects of SIP also has sufficient predictive skill in WRF.The main difference lies again at higher temperatures, where the RaFSIPv1 results in slightly overestimated LWC compared to the ALLSIP predictions, especially during fall and winter.

Cloud Radiative Forcing at the Surface
The preceding discussion illustrates how SIP can significantly alter the number concentration and size distribution of ice and liquid particles in polar MPCs.Such changes in the macro-and microphysical properties of clouds are expected to impact the two most critical cloud feedbacks associated with the polar regions: cloud optical depth and cloud-sea-ice feedbacks (e.g., Goosse et al., 2018).For instance, changes in the albedo of MPCs resulting from an increase in the amount of small liquid droplets contained in them can lead to enhanced shortwave radiation reflected back to space and reduced warming at the surface (Murray et al., 2021;Tan & Storelvmo, 2019).During the polar night, however, an increase in low-level cloud coverage has been found to increase the longwave radiation emitted to the surface that is trapped as heat due to the stable stratification conditions in the Arctic, ultimately resulting in surface warming (Ebell et al., 2020;Kay et al., 2016).In this section, we aim to evaluate the performance of the RaFSIP parameterization concerning the anomalies observed in the simulated cloud radiative forcing at the surface (CRF surf ) with respect to the ALLSIP simulation (ΔCRF surf , Figure 8).Calculations of CRF surf and ΔCRF surf are elaborated in supporting information of Young et al. (2019), which also accounts for additional modifications in the shortwave fluxes to account for the highly reflective Arctic surface (Vavrus, 2006).and locally up to ∼60 Wm 2 during fall and winter (Figures 8a and 8b).Ground-based remote sensing observations collected at Ny-Ålesund over two years (2016)(2017)(2018) found positive monthly net CRF surf between September and April/May, with a mean longwave warming of ∼50 Wm 2 (Ebell et al., 2020), which is in line with our results.In contrast, during summer, the higher occurrence of liquid in the simulated polar clouds (Figure 7n) combined with the lower surface albedo leads to a strong cloud-induced cooling at the surface by up to ∼ 100 Wm 2 (Figure 8d), again in agreement with Ebell et al. (2020).During spring, polar clouds have almost a neutral effect over the sea-ice and a slight cooling effect over the open ocean in the lower latitudes (Figure 8c).
The anomalies in longwave and shortwave radiation are calculated by the mean differences in the corresponding predictions between RaFSIP and ALLSIP, the combination of which produced the final surface anomalies, ΔCRF surf (Young et al., 2019).Note that the mean ΔCRF surf , is projected onto a 49 × 49 mesh by averaging over 3 × 3 sub-grids, similar to the averaged maps shown in Figure 4.During spring and summer, for instance, red colors in ΔCRF surf denote an increase in shortwave radiation reaching the ground in the RaFSIP simulation, implying the presence of more glaciated clouds than in ALLSIP, while blue colors signify a decrease in shortwave radiation suggesting fewer glaciated clouds than in ALLSIP.Replacing the detailed inline parameterizations of SIP in the M09 scheme with the new RaFSIP parameterization can lead to CRF surf biases between 2.5 and 2.5 Wm 2 , with a non-uniform pattern across the model domain (Figure 6e-6h).Compared to the CRF surf magnitudes illustrated in Figure 8, such radiative biases can be considered negligible.
During summer, when ice enhancement from SIP is the strongest, largest ΔCRF surf values are largely observed over the sea-ice and ocean grid cells at latitudes higher than 70°N (Figure 8h).At this time, RaFSIPv2 predictions of the ice (Figure 6p) and liquid (Figure 7p) cloud phase partitioning can be biased by up to a factor of 3, leading to an underestimation of the liquid-dominating clouds compared to ALLSIP, and consequently, an excess of up to ∼2.5 Wm 2 of shortwave radiation at the surface.However, in the continental and ocean regions of lower latitudes, mostly negative CRF surf biases are observed, indicating that RaFSIPv2 is lacking some ICNCs compared to ALLSIP there.A similar ΔCRF surf pattern characterizes spring (Figure 8g), with the colors shifted toward lighter tones, indicating relatively smaller radiative errors (between 1.5 and 1.5 Wm 2 ).Fall exhibits a cooling effect (down to 1 Wm 2 ), with some localized areas of warming (∼1 Wm 2 ) spread around the WRF domain (Figure 8e), while winter is the season with the minimum SIP contribution and the smallest cloud radiative biases, confined mostly to the bottom part of the model domain (Figure 8f).
Figure S15 in Supporting Information S1 displays the WRF outputs when RaFSIPv1 replaces RaFSIPv2 parameterization.The obtained ΔCRF surf ranges from 3 to 3 Wm 2 , indicating slightly higher radiative biases in all seasons compared to the results shown in Figure 8.The differences between the two SIP parameterizations are most significant in the warm season (Figures S15g and S15h in Supporting Information S1), where all positive biases resulting from RaFSIPv2 are replaced by negative ones.The most considerable negative bias appears in the sea area to the east of Greenland (latitudes between 70 and 80°N), where ∼3 Wm 2 less shortwave radiation reaches the surface during summer (Figure S15h in Supporting Information S1) compared to the ALLSIP simulation, indicating the tendency of RaFSIPv1 to overestimate the amount of liquid hydrometeors in summer MPCs.Its reduced predictive ability to capture all the enhanced SIP events present in the ALLSIP simulation stems mainly from the fact that it is limited by the presence of PIP, and hence it cannot be activated for the whole range of conditions that would favor ice multiplication.

Summary and Outlook
This study introduces a new ML-based framework aimed at representing the impact of ice multiplication in stratiform MPCs within large-scale models.The new RF parameterization, referred to as RaFSIP, is developed using two different approaches to describe the effect of SIP at temperatures as low as 25°C.The first approach indirectly describes SIP through the use of an ice multiplication factor applied to the primary ice production rates, which makes this version more linked with the PIP parameterizations employed in M09.The second approach directly predicts the SIP tendencies, providing a more explicit representation of the underlying physical processes.RaFSIP is trained on two years of 10-km horizontal grid spacing WRF simulations over the pan-Arctic region.We evaluate the parameterization both offline and online, and demonstrate that RaFSIP can accurately predict tendencies of unresolved physics in the 50-km horizontal grid spacing domain, comparable to the conventional SIP schemes.RaFSIP skillfully captures how SIP affects the glaciation state and evolution of polar stratiform MPCs, predicting their ice-and liquid-phase partitioning with an uncertainty of less than a factor of 3. Importantly, RaFSIP exhibits only minor radiative biases, not exceeding ±3 Wm 2 at the climatologically sensitive Arctic surface.This implies that biases in the predicted SIP rate or IEF (depending on the version of RaFSIP used) do not translate into significant errors in simulated cloud microphysical properties and cloud radiative forcing.
The intrinsic property of RF algorithms to make conservative predictions by averaging over the predictions of the individual decision trees allowed the coupled WRF-RaFSIP scheme to run efficiently and without instabilities for a year.Unlike other ML algorithms designed primarily to speed up computationally expensive microphysical processes in models (e.g., Gettelman et al., 2021), WRF-RaFSIP demonstrates similar computational costs to simulations involving detailed SIP microphysics.The primary advantage of our approach lies in offering a streamlined and easily implementable solution for representing SIP in large-scale models, overcoming challenges associated with implementing complex microphysics at coarse spatial scales.This is particularly significant since important SIP processes are still absent from the default releases of NWP and climate models.Within RaFSIP, the effect of the three most prolific SIP mechanisms is described using inputs, such as temperature, ice-and liquid water content, riming tendencies, and relative humidity with respect to ice, which are typically present in most state-of-the-art models.This can greatly facilitate its integration in host models with horizontal grid spacings comparable to or larger than the 10-km grid size used for training RaFSIP.Given the established impact of SIP on cloud evolution, rapidly glaciating them, altering precipitation and eventually climate (Atlas et al., 2022;Sotiropoulou et al., 2024;Zhao & Liu, 2021), the ability to incorporate its representation in large-scale models through a simplified approach is crucial.This marks a significant step toward developing advanced ML-based approaches for cloud microphysics that are either scale-adaptive or more computationally efficient than traditional physical parameterizations.
The incorporation of the RaFSIP framework in the microphysics scheme of 3 GCMs, namely ECHAM-HAM, NorESM2, and EC-Earth, participating in the CMIP assessments is one of the undergoing studies we are currently addressing (Costa-Surós et al., 2023;Ickes et al., 2023).The modular nature of the new SIP parameterization allows for seamless integration with the dynamical core of each GCM model, provided that the input features are well located inside the stratiform cloud microphysics routine.Preliminary results indicate that utilizing grid-box averaged atmospheric variables predicted by GCMs as inputs to the RaFSIP scheme leads to computationally efficient simulations without numerical instabilities or explosive SIP (not shown), further demonstrating that RaFSIP can work stably across scales.This model intercomparison project will evaluate the robustness of our approach, pointing out any potential limitations.
As the RaFSIP scheme is trained on a wide range of stratiform conditions where SIP is likely to occur, it can effectively capture the underlying mechanisms behind ice multiplication in any type of single or multi-layered clouds in the polar regions.The limitations for our work stem mostly from the fact that RaFSIP is trained based on regional simulations, and it therefore may not be applicable for the entire globe (e.g., tropical marine or deep convective clouds).For convective clouds, a notable limitation arises from the absence of vertical velocity as an input parameter to RaFSIP-crucial for determining the duration cloud particles spend within the SIPfavorable temperature regime.Recognizing the importance of cloud dynamics for SIP, training RaFSIP using outputs from higher-resolution km-scale grids or even Large Eddy Simulation models would be imperative.In the context of growing convective clouds, where the HM process may surpass BR in efficiency (Waman et al., 2022), and DS may be more active due to the presence of larger raindrops, potential adaptations for a scale-adaptive RaFSIP scheme (Chen et al., 2023) could be more urgent than in stratiform clouds.This is because of the scale dependence of LWC, which can be an important input for predicting the effect of HM.Additionally, in fast convective downdrafts, fragmentation during ice particle sublimation under ice subsaturated cloudy conditions (Deshmukh et al., 2022) might also be an important source of SIP particles following BR, as reported by Waman et al. (2022).This implies that additional SIP processes should be included in the RaFSIP training data set to extend its applicability in different regions worldwide, where the importance of SIP is recognized (Hoose, 2022).
Another limitation of the RaFSIP scheme arises from its training on a specific bulk microphysics scheme within WRF (M09).Particularly for the RaFSIPv1 scheme, the dependence of the IEF on the specific temperaturedependent primary ice nucleation parameterizations of M09 raises concerns about its applicability in other microphysical schemes.To address this limitation, the more versatile RaFSIPv2 scheme was developed.By utilizing total IWC as an input variable, constrained by temperature and RHI, RaFSIPv2 can more effectively capture the favorable SIP conditions independent of the employed PIP parameterizations.Additionally, the detailed representation of BR in M09 required adopting prescribed rimed fraction and ice habit in the 2-year Journal of Advances in Modeling Earth Systems 10.1029/2023MS003923 training data set, which may introduce uncertainties in the predicted secondary ice fragments.Using a more sophisticated microphysics scheme in WRF that predicts rime mass fraction (e.g., the "P3"scheme, Milbrandt et al., 2021;Morrison & Milbrandt, 2015) could mitigate these uncertainties.The absence of explicit ice particle habit description (e.g., Jensen et al., 2017) is not expected to significantly impact the representation of BR in Arctic clouds, as shown by (Sotiropoulou, Ickes, et al., 2021).Despite these limitations, our analysis demonstrates that, by leveraging a relatively small number of grid points from a regional climate modeling framework, the training data set of RaFSIP can be readily expanded to address inherent caveats in our approach.
Beyond its modeling applications, the adaptable nature of the RaFSIP framework holds potential for diagnosing SIP contributions from field observations.This applies specifically to the RaFSIPv1 scheme, as the IEF can be derived from either in-situ or ground-based remote sensing observations (Wieder et al., 2022), provided that simultaneous measurements or retrievals of ICNCs and INPs are available.An observational challenge in quantifying IEF due to SIP is ensuring that MPCs remain unaffected by snow particles from the surface or seeding from higher-level clouds.Our study has already outlined the essential input variables for constraining IEF.Such an observation-constrained version of RaFSIP could offer data sets for model evaluation, providing valuable insights into SIP across various regions and cloud types.
This study serves as a proof of concept for the potential of ML techniques and fine-resolution data for the development of computationally efficient and accurate subgrid-scale parameterizations of cloud microphysical processes, addressing the current gaps and limitations in large-scale models.By integrating physics-guided ML algorithms and utilizing feature importance metrics, we increase the transparency of the inner workings of these algorithms, even if they are initially considered uninterpretable (Beucler et al., 2021;McGovern et al., 2019).This is a promising pathway toward model simplification and improved process representation within GCMs.Moving forward, further exploration and implementation of such techniques can contribute to advancements in GCMs, ultimately leading to more robust and comprehensive climate simulations.
We assume that collisions between ice particles not resulting in the collection of the smaller particle by the more massive one can lead to ice fragmentation after BR.The efficiency of the BR process is therefore assumed to be equal to 1-E col .
To parameterize the number of fragments generated after ice-ice collisions in M09, we followed the formulation of Phillips et al. (2017), developed using the experimental observations by Vardiman (1978) and Takahashi et al. (1995), further considering the physics of collisions.Among the possible collisions between ice particles in M09, BR is allowed only in cases where the smallest particle (or the less dense one) has some rime on its surface from accreted cloud droplets or raindrops.This step is crucial for avoiding overestimation of the number of fragments generated, given that both experimental studies involved large, rimed ice particles (Georgakaki et al., 2022).The formula describing the number of secondary ice fragments per collision is the following: where: ) 2 is the initial kinetic energy (in J), with m 1 and m 2 being the masses of the colliding particles (in kg), while |Δu n12 | = {(1.7 u n1 u n2 ) 2 + 0.3 u n1 u n2 } 1/ 2 is a correction term proposed by Mizuno et al. (1990) and Reisner et al. (1998), to account for underestimates when the colliding ice particles have similar terminal velocities (i.e., u n1 ≈ u n2 ) • a is the surface area of the particle undergoing fragmentation (in m 2 ), defined as a = πD 2 , with D its size (in m) • A represents the number density of breakable asperities on the colliding surfaces (in m 2 ) • C is the asperity-fragility coefficient (in J 1 ), which is empirically derived to account for different collision types • γ is a parameter related to the riming intensity The above mentioned parameters are adapted based on the type of ice collisions.For collisions involving break-up of planar cloud ice or snow particles.
• A = a 0 3 + max( 2a 0 3 a 0 9 |T 258|,0) , in which a 0 = 3.78 × 10 4 (1 + 0.0079 D 1.5 ) • γ = 0.3 • C = 6.3 × 10 6 ψ Multiplying the number of SIP fragments (N BR ) with the number tendency of collisions between the corresponding ice particles yields the production rate of secondary ice splinters due to BR, dNi BR dt (in kg 1 s 1 ) that will be added to the conservation equation of cloud ice particles.Note that, an upper limit of N BR equal to 100 fragments is also imposed, as suggested by Phillips et al. (2017).

A3. Droplet-Shattering
To incorporate the DS mechanism in M09 we adopted the parameterization developed by Phillips et al. (2018), which considers two different modes, each represented by a different mathematical formulation, describing the freezing of supercooled raindrops.In "mode 1," the freezing of the raindrop is initiated after a collision with a smaller ice particle or an INP (immersion freezing mode), while "mode 2" requires a collision with a more massive ice particle, such as snow or graupel.
In mode 1, the freezing raindrop maintains its spherical symmetry, and upon freezing, it can generate both tiny and big fragments, with its efficiency peaking around 15°C.The formulation for mode 1 is derived by pooling multiple laboratory data sets into a Lorentzian function of temperature and a polynomial expression of the drop size.The mathematical formulation describing the total number of fragments (N) is as follows: where.
• T is the temperature (in K) • D r is the size of the freezing raindrop (in mm) • Ξ(D r ) and Ω(T ) are cubic interpolation functions, preventing the activation of DS for D r < 0.05 mm and T > 3°C The number of large fragments (N B ) is: The parameters ζ, η, Τ 0 , β, ζ Β , η Β , Τ Β,0 , are described in Phillips et al. (2018).It is only the tiny splinters, calculated as N T = N-N B , that are added to the cloud ice category.The larger fragments are initiated in the model as graupel, snow, or frozen drops, depending on the mechanism that triggered their formation.
In mode 2, the collision with more massive ice particles disrupts the spherical symmetry of the raindrop, generating only tiny fragments upon freezing.Although there is only one dedicated laboratory study on this freezing mode (James et al., 2021), Phillips et al. (2018) proposed a theoretical, energy-based formulation to represent the number of tiny splinters per drop accreted: Here, DE is the dimensionless energy, given as the ratio of the initial kinetic energy (K 0 ) to the surface energy.DE crit is the critical value for the onset of splashing upon impact.The parameter f(T) represents the initial frozen fraction of a supercooled drop, and Φ(T) is an empirical fraction representing the probability of a secondary drop containing a frost secondary ice particle.Further details about the empirical parameters and uncertainties underlying the mathematical formulations are discussed in Phillips et al. (2018).

Figure 1 .
Figure 1.Schematic diagram of the methodology followed to develop the two versions of the RaFSIP parameterization.The map on the left shows the two weather research forecasting (WRF) domains used to simulate the 2-year period required for training RaFSIP.The 10-km horizontal grid spacing domain is outlined by purple dashed lines, while colors indicate terrain heights (left axis) and sea ice concentration (right axis).The 4 triangles represent the regions from where the WRF data set is extracted.This 2-year data set is split into ∼83% training data set (20 months) and ∼17% testing data set (4 months).The RaFSIPv1 (RaFSIPv2) parameterization is trained to predict the ice enhancement factor (SIP rate ) as well as the transfer of mass from cloud droplets (Qc tr ) or raindrops (Qr tr ) to cloud ice due to SIP.The snapshot corresponds to the monthly average of January 2016.

Figure 2 .
Figure 2. Schematic illustrating the sequential steps within the M09 scheme of the Weather Research Forecasting model, leading to the activation of different random forest regressor (RFR) models, represented by the dark blue boxes in the final stage of the flowchart.Each RFR model is selectively invoked to predict the impact of various secondary ice production (SIP) mechanisms based on prevailing meteorological and microphysical conditions.The inset, delineated by black dashed lines, specifically pertains to the RaFSIPv2 parameterization, where the effect of SIP can also be considered in the warm regime between 0 and 3°C.

Figure 3 .
Figure 3. Normalized histograms (i.e., frequency is scaled by the total number of predictions) comparing true Weather Research and Forecasting (WRF) model outputs with forestALL predictions in offline tests focusing on the 10-km horizontal grid spacing WRF data.Top panel: RaFSIPv1 predictions of the ice enhancement factor due to (a) BR, (b) HM, (c) DS, and the transferred masses from (d) cloud droplets (Qc tr ) and (e) raindrops (Qr tr ) to cloud ice.Bottom panel: RaFSIPv2 predictions of the secondary ice production rate due to (f) BR, (g) HM, (h) DS, and the transferred masses (i) Qc tr and (j) Qr tr .The black line represents the one-to-one line.

Figure 4 .
Figure 4. Latitude-longitude maps display median ice crystal number concentration from the 50-km horizontal grid spacing domain of the Weather Research and Forecasting model from the (a, e, i, m) CONTROL, (b, f, j, n) ALLSIP, and (c, g, k, o) RaFSIP simulations.Corresponding (d, h, l, p) R 2 scores between ALLSIP and RaFSIPv2 simulation outputs are also shown, with R 2 maps spatially averaged over 3 × 3 grid cell regions.Statistics are calculated for the 4 seasons: fall (top panel), winter (second panel), spring (third panel) and summer (bottom panel).

Figure 5 .
Figure 5. Median ice crystal number concentration values presented for distinct latitudinal zones, temperature ranges, and seasons, extracted from the 50-km horizontal grid spacing domain of the Weather Research and Forecasting model from the CONTROL (black), ALLSIP (green), and RaFSIP (purple) simulations.Cold period values (a, c, e) span from September 2019 to February 2020, while warm period values (b, d, f) cover March-August 2020.Markers distinguish whether the median ice crystal number concentration (ICNC) is derived from the entire latitude zone (open circle markers) or from glaciated clouds (glaciation fraction >0.98) over sea-ice regions (sea-ice concentration >50%) (filled rhombus markers).The gray rhombuses represent the median ICNC derived from the 10-year satellite data set presented in Figure 1 of PP22.Error bars represent the 95% confidence intervals for glaciated clouds over sea-ice.

Figure 6 .
Figure 6.Bivariate joint PDF of median ice crystal number concentration (ICNC) defined in terms of both temperature and glaciation fraction from the 50-km horizontal grid spacing domain of the Weather Research and Forecasting model from the CONTROL (a, e, i, m), ALLSIP (b, f, j, n) and RaFSIPv2 (c, g, k, o) simulations.The normalized histograms of the 2D-binned median ICNC outputs from the ALLSIP simulation are also plotted versus the RaFSIPv2 predictions (d, h, l, p).The black line represents the one-to-one line, while the gray dashed lines delimit the area where the ALLSIP values are over or under-estimated by a factor of 3. Statistics are calculated for the 4 simulated seasons: fall (top panel), winter (second panel), spring (third panel) and summer (bottom panel).

Figure 7 .
Figure 7. Bivariate joint PDF of median liquid water content (LWC) defined in terms of both temperature and glaciation fraction from the 50-km horizontal grid spacing domain of the Weather Research and Forecasting model from the CONTROL (a, e, i, m), ALLSIP (b, f, j, n) and RaFSIPv2 (c, g, k, o) simulations.The normalized histograms of the 2D-binned median LWC outputs from the ALLSIP simulation are also plotted versus the RaFSIPv2 predictions (d, h, l, p).The black line represents the one-to-one line, while the gray dashed lines delimit the area where the ALLSIP values are over or under-estimated by a factor of 3. Statistics are calculated for the 4 simulated seasons: fall (top panel), winter (second panel), spring (third panel) and summer (bottom panel).

Figure 8
Figure8illustrates the annual cycle of the cloud impact on the surface radiative energy budget, with blue and red shades indicating cloud-induced cooling and warming, respectively.Polar clouds have a strong longwave greenhouse warming effect during the non-summer months, with maximum CRF surf reaching up to ∼40 Wm 2

Figure 8 .
Figure 8. Latitude-longitude maps of (top panel): seasonally averaged cloud radiative forcing at the surface (CRF surf ) derived from the 50-km horizontal grid spacing domain of the Weather Research and Forecasting (WRF) model from the ALLSIP simulation, and (bottom panel): anomalies with respect to the ALLSIP simulation (ΔCRE surf ) derived from the WRF simulation coupled with RaFSIPv2, extracted respectively for (a, e) fall, (b, f) winter, (c, g) spring and (d, h) summer.Note that ΔCRE surf is spatially averaged over 3 × 3 grid cell regions.

Table 1
Root-Mean-Square Error (RMSE) of RaFSIPv1 Parameterization Predictions When Compared Against True WRF Outputs in Offline Tests Note.The corresponding coefficient of determination (R 2 ) values are also shown in parentheses.The offline analysis focuses on the 10-km horizontal grid spacing data, while scores are based on normalized outputs.

Table 2
RMSE of RaFSIPv2 Parameterization Predictions When Compared Against True WRF Outputs in Offline TestsNote.The corresponding R 2 values are also shown in parentheses.The offline analysis focuses on the 10-km horizontal grid spacing data, while scores are based on normalized outputs.