In the current article, the problem of in-flight ice accumulation on multi-element airfoils is studied numerically. The analysis starts with flow field computation using the Hess-Smith panel method. The second step is the calculation of droplet trajectories and droplet collection efficiencies. In the next step, convective heat transfer coefficient distributions around the airfoil elements are calculated using the Integral Boundary-Layer Method. The formulation accounts for the surface roughness due to ice accretion. The fourth step consists of establishing the thermodynamic balance and computing ice accretion rates using the Extended Messinger Model. At low temperatures and low liquid water contents, rime ice occurs for which the ice shape is determined by a simple mass balance. At warmer temperatures and high liquid water contents, glaze ice forms for which the energy and mass conservation equations are combined to yield a single first order ordinary differential equation, solved numerically. Predicted ice shapes are compared with experimental shapes reported in the literature and good agreement is observed both for rime and glaze ice. Ice shapes and masses are also computed for realistic flight scenarios. The results indicate that the smaller elements in multielement configurations accumulate comparable and often greater amount of ice compared to larger elements. The results also indicate that the multi-layer approach yields more accurate results compared to the one-layer approach, especially for glaze ice conditions.