Accurate simulation of icing is important for the assessment of several potential icing scenarios and complex icing regulations. However, performing all possible icing scenarios is a demanding process in terms of computational cost, especially when modification of the geometry due to ice accretion is required. Additionally, aircraft icing safety assessment necessitates an evaluation of the accumulated ice. Thus, numerical representation of the non-linear and complex geometries is essential for the parametrization of this ice. Indeed, surrogate models have the capability of predicting these complex, non-linear shapes. For this purpose, a method for ice accretion prediction on a selected airfoil, NACA 22112, is proposed in this study with different surrogate models that will later be used for fast prediction in 6DOF simulations to directly evaluate its effects on aerodynamic performance during flight. The required datasets in order to train for clean and iced airfoils are based on numerical analysis results obtained through the FENSAP-ICE 2022 R1 commercial tool with a multi-shot technique. They are generated by varying four variables (liquid water content, ambient temperature, median volumetric diameter, and exposure time), which are the most prominent atmospheric or cloud parameters for ice shapes. The combination of these input datasets is selected based on the 14 CFR Part 25 Appendix-C envelopes, and ice shapes are modeled by applying the Fourier series expansion approach. According to the results, nearly 30 Fourier coefficients can accurately capture nonlinear rime ice shapes within acceptable deviations. Moreover, surrogate models such as artificial neural networks and Gaussian processes are compared to predict these coefficients in terms of their ability to capture targeted ice shapes.