Ice accretion threatens aircraft safety. With the wide application of unmanned aerial vehicles (UAVs), the design of ice-tolerant UAVs has become a problem that must be solved. Forty conditions for the continuous/intermittent maximum icing conditions were sampled in Appendix C of Federal Aviation Regulations Part 25. Herein, numerical simulations of icing were performed on an NACA 0009 airfoil for 5, 15, and 30 min, and the ice mass and ice shapes were obtained at different times. Numerical simulations of the aerodynamic characteristics of the iced configuration at 5, 15, and 30 min were conducted, and the coefficients of lift, drag, and pitch moment were obtained. Surrogate models of the ice shape, mass of the ice accretion, lift coefficient, drag coefficient, and pitch moment coefficient at different moments were built based on proper orthogonal decomposition and kriging interpolation. The results demonstrate that the surrogate models accurately predicted the ice shape, ice mass, lift coefficient, drag coefficient, and pitch moment coefficient at different moments. Compared with the numerical simulation results, the maximum relative errors of the ice mass, lift coefficient, drag, and pitch moment predicted by the surrogate models were 7.8%, 3.4%, 3.9%, and 7.6%, respectively. This method can help in designing the ice-tolerant UAVs and envelope determination under icing conditions.