The design of efficient thermal ice protection systems is a challenging task as these systems operate in complex environments involving several coupled physical phenomena such as phase change, boundary-layer flow, and heat transfer. Moreover, certification rules are becoming more stringent, and there is a strong incentive for the reduction of fuel consumption. In this context, numerical tools provide a powerful asset during the design phase but also to gain insight into the physical mechanisms at play. This article presents modeling and simulation strategies for thermal ice protection systems. First, the model describing the behavior of the thermal protection system is presented. Second, a model and associated numerical method is presented for unsteady ice accretion. Third, the coupling methodology between the ice accretion solver and the heat conduction solver is described. In the fourth part, different methods to simulate the boundary-layer flow are described. Finally, some relevant examples are presented, both in steady and unsteady configurations.