Numerically predicted roughness distributions obtained in in-flight icing simulations with a beading model are used in a quasi-steady manner to compute ice shapes. This approach, called "Multishot," uses a number of steady flow and droplet solutions for computing short intervals (shots) of the total ice accretion time. The iced geometry, the grid, and the surface roughness distribution are updated after each shot, producing a better match with the unsteady ice accretion phenomena. Comparisons to multishot results with uniform roughness show that the evolution of the surface roughness distribution has a strong influence on the final ice shape. The ice horns that form are longer and thinner compared to constant roughness results. The constant roughness approach usually fails to capture the formation of the pressure side horns and under-predicts the thickness of the ice in this region. With the variable roughness produced by the beading model, the ice shape in general is captured more accurately, better matching the experimental ice shapes. Results are provided for NACA 0012 and SC 1095 airfoils for 2D icing. The flow conditions range between Mach numbers of 0.3 to 0.6, free stream static temperatures of -27°C to -6°C, free stream liquid water contents of 0.5 g/m₃ to 1.4 g/m₃, Reynolds numbers of 2.6 to 4.4 million, and droplet diameter of 20 microns. Total ice accretion times range between 45 seconds to 7 minutes. The FENSAP-ICE icing simulation system is utilized to carry out the CFD-Icing calculations in 2D and 3D.