Computational Ice Dynamics represents one of the most challenging branches of fluid dynamics. Interest in this numerical discipline is driven by aerodynamic degradation from in-flight ice-accretion, a major flight safety concern. The development and assessment of advanced numerical methods have become essential to predict and mitigate these phenomena, reducing reliance on costly experimental campaigns in large ice wind tunnels (IWT) [1]. Among emerging techniques, immersed boundary (IB) methods coupled with Level-Set formulations represent promising solutions to overcome mesh-generation bottlenecks and enable high-fidelity simulations within industrial time constraints. Compared to Lagrangian approaches, which explicitly track the ice-front and are favored for their simplicity, the Level-Set Method (LSM) offers significant advantages in accuracy and robustness especially for multi-step analyses [2]. Lagrangian methods often suffer from mass conservation errors and geometric inconsistencies, particularly in concave regions, due to self-intersecting meshes and complex topology changes. In contrast, LSM represents the interface implicitly on a fixed Eulerian grid, automatically handling ice-front formation without requiring costly mesh re-generation at each time-step. The present research concerns a multi-step process that couples a Cartesian IB-method (solving air and water phases) with a Eulerian ice-front tracking based on a level-set function. The latter consists of a signed distance from the wall which is computed in a pre-processing phase by solving the Eikonal equation. At each time-step, the thermodynamic water-film solver provides the so-called “ice-front velocity” on the wall surface whose module is linked with the ice-mass growth rate. A conservative PDE is applied to propagate this velocity in the domain volume to allow the advection of the level-set function. Periodic reinitialization of the level-set field ensures numerical stability and accuracy. Overall, the proposed method combines geometric flexibility with computational efficiency, making it ideal for aero-icing simulations where interface evolution and topology changes are critical.