Representing single or multi-layered mixed-phase clouds (MPCs) accurately in global climate models (GCMs) is critical for capturing climate sensitivity and Arctic amplification. Ice multiplication, or secondary ice production (SIP), can increase the ice crystal number concentration (ICNC) in MPCs by several orders of magnitude, affecting cloud properties and processes. Here, we propose a machine-learning approach, called Random Forest SIP (RaFSIP), to parameterize the effect of SIP on stratiform MPCs. The RaFSIP scheme uses few input variables available in models and considers rime splintering, ice-ice collisional break-up, and droplet-shattering, operating at temperatures between 0 and -25 ˚C. The training dataset for RaFSIP was derived from two-year pan-Arctic simulations with the Weather Research and Forecasting (WRF) model with explicit representations of SIP processes. The RaFSIP scheme was evaluated offline against WRF simulation outputs, then integrated within WRF. The parameterization exhibits stable performance over a simulation year, and reproduced predictions of ICNC with explicit microphysics to within a factor of 3. The coupled WRF-RaFSIP scheme can replicate regions of enhanced SIP and accurately map ICNCs and liquid water content, particularly at temperatures above -10 ˚C. Uncertainties related to the RaFSIP representation of MPCs marginally affected surface cloud radiative forcing in the Arctic, with radiative biases of lower than 3 Wm-2 compared to simulations with explicit SIP microphysics. Training from a few high-resolution model grid points did not limit the predictive skill of RaFSIP, with the approach opening up new avenues for model simplification and process description in GCMs by physics-guided machine learning algorithms.