Steady and unsteady three-dimensional viscous turbulent aero-icing simulations are computationally expensive, especially for certification campaigns when broad parametric studies are needed. To overcome the computational effort of such investigations, a Reduced Order Modeling approach, based on Proper Orthogonal Decomposition and Kriging interpolation, is proposed. Using a database of high-fidelity numerical simulations, experimental data, or combinations of both, the proposed technique allows approximating solutions by linear combination of a limited number of eigenfunctions. Bayesian Kriging, a recent variant, is used to obtain the scalar coefficients for the expansion. The accuracy of the proposed method is assessed against reference solutions from two- and three-dimensional aero-icing simulations as well as against experimental data.