The in-flight ice accretion simulations are typically performed using a quasi-steady formulation through a multi-step approach. As the ice grows, the geometry changes, and an adaptation of the fluid volume mesh used by the airflow and droplet-trajectory solver is required. Re-meshing or mesh deformation are generally employed to do that. The geometries formed are often complex ice shapes increasing the difficulty of the re-meshing process, especially in three-dimensional simulations. Consequently, difficulties are encountered when trying to automate the process. Contrary to the usual body-fitted mesh approach, the use of immersed boundary methods (IBMs) allows solving, or greatly reducing, this problem by removing the mesh update, facilitating the global automation of the simulation. In the following paper, an approach to perform the airflow and droplet trajectory calculations for three-dimensional simulations is presented. This framework utilizes only immersed boundary methods. In particular, two methods are presented. On the one hand, a ghost-cell Immersed Boundary approach has been developed to solve the aerodynamics. The Euler equations are solved at fluid points, whereas the solution is forced in the vicinity of the obstacle at some particular cells (IB target points), in order to mimic a slip boundary condition. Special attention has been given to the applied boundary conditions as well as to the location of these IB target points. In fact, instead of the most commonly used approaches where the IB target points are placed in the solid or in the fluid, the case where these IB target points lie astride the obstacle (namely “GC Surrounding” in the following), in the fluid and solid regions, is studied. On the other hand, to solve the droplet trajectory equations, the penalization method, already present in IGLOO2D (the 2D ice accretion suite developed at ONERA), has been expanded to the three-dimensional simulations. The two immersed boundary methods are compared and numerically tested on different cases. First, a mesh refinement study is performed for weakly compressible flow around a cylinder. In this case, the solution is compared with that obtained using a body-fitted simulation and it serves as a verification of the method. Next, the approach is used in two other different cases. The first case involves an iced GLC305 airfoil, which is characterized by its complex geometry, under compressible subsonic conditions. The second case is a three-dimensional simulation in which the presented approach is used to analyze the weakly compressible flow around a sphere.