In-flight ice accumulation on aircraft components poses significant safety concerns and must be care-fully considered during design cycles. Aircraft exposed to clouds of supercooled water droplets under specific atmospheric conditions are susceptible to ice formation. Current icing prediction tools struggle to accurately reproduce the complex ice morphologies observed experimentally, such as ice feathers or scallops. Since the aerodynamics, droplet field, and thermodynamic balance depend on the ice shape, ice accretion must be simulated in steps. Consequently, robust mesh generation tools are required, and mesh complexity remains a major limitation. To address this issue, Immersed Boundary Methods (IBM) have been been explored for icing phenomena, notably in previous work by Capizzano [1]. The proposed method implements IBM to compute both the flow field and the droplet field, while accounting for the thermodynamic balance at the boundary. The flow field is solved using a ghost-cell approach [2], whereas the droplet field is computed in an Eulerian framework with a penalization method [3]. Wall functions are incorporated to enable simulations at high Reynolds numbers. This approach improves mesh robustness in the IBM solver by relying on an automatically generated mesh, locally refined near the ice geometry through heuristic procedures. Compatible libraries such as CGAL [4] are integrated into the in-house software Chapel Multi-Physics Simulation (CHAMPS), developed at Polytechnique Montr´eal. Ice fronts obtained from a stochastic icing model [5] are re-injected into the IBM Reynolds-averaged Navier-Stokes (RANS) solver, and ice accretion is generated layer by layer. Results will show multi-layer ice shapes on industrial configurations. To the authors’ knowledge, this work presents the first model coupling IBM with a multi-layer stochastic ice accretion framework, contributing to improved numerical icing prediction for aircraft geometries.