Review of Wind Turbine Icing Modelling Approaches

: When operating in cold climates, wind turbines are vulnerable to ice accretion. The main impact of icing on wind turbines is the power losses due to geometric deformation of the iced airfoils of the blades. Signiﬁcant energy losses during the wind farm lifetime must be estimated and mitigated. Finding solutions for icing calls on several areas of knowledge. Modelling and simulation as an alternative to experimental tests are primary techniques used to account for ice accretion because of their low cost and effectiveness. Several studies have been conducted to replicate ice growth on wind turbine blades using Computational Fluid Dynamics (CFD) during the last decade. While inﬂight icing research is well developed and well documented, wind turbine icing is still in development and has its peculiarities. This paper surveys and discusses the models, approaches and methods used in ice accretion modelling in view of their application in wind energy while summarizing the recent research ﬁndings in Surface Roughness modelling and Droplets Trajectory modelling. An An additional section discusses research on the modelling of electro-thermal icing protection systems. This paper aims to guide researchers in wind engineering to the appropriate approaches, references and tools needed to conduct reliable icing modelling for wind turbines.


Introduction
Due to the rising concern over climate change, actual global trends move increasingly toward reducing greenhouse gas emissions. Among these actions is the investment in the exploitation of clean energy resources. According to Natural Resources Canada, about 78% of global greenhouse gas emissions from human activity are from production and energy consumption [1]. Consequently, during the last few decades, investment in clean energy technologies has been growing progressively. Wind energy production has been booming to become a significant part of the energy market and research worldwide [2]. The "BP Statistical Review of World Energy2021" revealed that wind provided the most significant contribution to renewables' global growth in 2020 [3]. In Canada, electricity generation from renewable sources increased by 16% between 2010 and 2018, and wind energy has the most considerable growth [1]. Kumar,et al. [4] discussed the trends and the available technologies developed for wind energy in conjunction with their applications and devices of operation [4]. However, this technology is highly dependent on its geographic location.
In the Nordic countries, the wind potential is very high in winter; the wind is stronger, and air density is higher. Therefore, countries with a high wind potential and intensive cold climates aim to promote the exploitation of wind power resources [5]. The cumulative wind capacity installed in northern climates throughout Scandinavia, North America, Europe and Asia was around 127 GW at the end of 2015, with an annual predictable growth rate of 11.7 GW by 2020 [6]. However, for operating in severely cold climates, wind technology must adapt to challenging conditions with frequent low-temperature events and, in particular, icing. Figure 1 shows wind turbines in icing conditions in Canada. In the photo appears the largest vertical axis wind turbine. @Fahed Martini. Wind Turbine Icing is the phenomenon of ice build-up mainly on wind turbine blades. It is a common problem in Nordic countries, which has severe consequences on the operation and the life cycle of wind turbines. When wind turbine blades are subjected to supercooled water droplets under specific meteorological conditions, ice starts to build forward in different forms, types, and severities. The deposition of ice on the blades may affect the turbine's structural integrity, causing much damage to the turbines and the surrounding areas. Safety hazards have been reported, such as the throwing of large chunks of ice may threaten the safety of the public, employees, and nearby roads and facilities. [7]. Various types of icing consequences have been stated by the cited sources [8][9][10][11][12]. However, it is the shape rather than the mass or dimensions of ice that significantly impact wind energy production [7]. The accumulation of ice on the leading edge of blade section airfoils, even in small quantities, results in significant consequences on power production of the wind turbine due to the degradation of the aerodynamics performance associated with the iced airfoil deformation [7].
Research into wind turbine icing is a developing discipline with limited available resources for investigations [13]. Given that the phenomenon of ice accumulation is physically complex, the optimization of wind turbines operating under icing conditions depends essentially on experiments. However, experimental tests and wind farm observations being relatively costly and having many operational difficulties; numerical simulation has increasing importance as a means of compliance. By reproducing the ice accretion progression on wind turbines blades, a numerical simulation approach makes it possible to conveniently provide information on the aerodynamic and energy production losses, energy required for de-icing for different wind turbine configurations and weather conditions [14,15]. As computational cost drops, the use of computational fluid dynamics In the photo appears the largest vertical axis wind turbine. @Fahed Martini. Wind Turbine Icing is the phenomenon of ice build-up mainly on wind turbine blades. It is a common problem in Nordic countries, which has severe consequences on the operation and the life cycle of wind turbines. When wind turbine blades are subjected to supercooled water droplets under specific meteorological conditions, ice starts to build forward in different forms, types, and severities. The deposition of ice on the blades may affect the turbine's structural integrity, causing much damage to the turbines and the surrounding areas. Safety hazards have been reported, such as the throwing of large chunks of ice may threaten the safety of the public, employees, and nearby roads and facilities [7]. Various types of icing consequences have been stated by the cited sources [8][9][10][11][12]. However, it is the shape rather than the mass or dimensions of ice that significantly impact wind energy production [7]. The accumulation of ice on the leading edge of blade section airfoils, even in small quantities, results in significant consequences on power production of the wind turbine due to the degradation of the aerodynamics performance associated with the iced airfoil deformation [7].
Research into wind turbine icing is a developing discipline with limited available resources for investigations [13]. Given that the phenomenon of ice accumulation is physically complex, the optimization of wind turbines operating under icing conditions depends essentially on experiments. However, experimental tests and wind farm observations being relatively costly and having many operational difficulties; numerical simulation has increasing importance as a means of compliance. By reproducing the ice accretion progression on wind turbines blades, a numerical simulation approach makes it possible to conveniently provide information on the aerodynamic and energy production losses, energy required for de-icing for different wind turbine configurations and weather conditions [14,15]. As computational cost drops, the use of computational fluid dynamics (CFD) for wind turbine design and analysis is becoming increasingly widespread. It results in a better understanding of the phenomenon and its consequences on wind turbines [16]. Thanks to the development of advanced icing simulation tools for aeronautics, the bibliographic study showed increasing interest in modelling and simulation research adapting these tools for wind turbine icing simulation during the last decade. The most commonly used icing simulation tools are LEWICE, developed by NASA and ANSYS FENSAP-ICE, developed by NTI Inc. The latter is often used recently in the literature for wind turbine icing due to its availability and integration with Ansys FLUENT and CFX. (see Figure 2). Based on what has been mentioned above, it is appropriate to review wind turbine icing modelling and provide engineers, analysts, or similar professionals in computational wind energy with icing simulation tools primarily validated with the operational conditions of wind turbines.
(CFD) for wind turbine design and analysis is becoming increasingly widespread. It results in a better understanding of the phenomenon and its consequences on wind turbines [16]. Thanks to the development of advanced icing simulation tools for aeronautics, the bibliographic study showed increasing interest in modelling and simulation research adapting these tools for wind turbine icing simulation during the last decade. The most commonly used icing simulation tools are LEWICE, developed by NASA and ANSYS FENSAP-ICE, developed by NTI Inc. The latter is often used recently in the literature for wind turbine icing due to its availability and integration with Ansys FLUENT and CFX. (see Figure 2). Based on what has been mentioned above, it is appropriate to review wind turbine icing modelling and provide engineers, analysts, or similar professionals in computational wind energy with icing simulation tools primarily validated with the operational conditions of wind turbines. Figure 2. A shot of a dynamic ice accretion simulation at the surface of a wing using ANSYS-FESAP-ICE. Image was extracted from ANSYS webinar titled "Accelerating the Aircraft Icing Certification Process Using Simulation".
Unlike the previous icing studies, this paper, in addition to highlighting the problem of wind turbine icing in northern climates, focuses on the actual state-of-the-art of ice modelling and simulation. It provides guidelines to the knowledge needed for modelling the phenomenon and examines the numerical methods used to analyze ice. It provides the best practices and recommendations for the simulation of ice accretion on wind turbines. After reviewing hundreds of references available in the literature on wind turbine icing, investigations have been carried out on the art of modelling ice accretion, from the applicable theory and various mathematical and physical models to the methodological approaches adopted to conduct reliable modelling in this domain. The objective is to present and recommend trustworthy modelling tools that are capable, via numerical simulations, of predicting and quantifying icing and eventually adapting the wind turbines to handle severe winter conditions.
The present paper is structured as follows: Section 2 summarizes the computational approaches employed in the simulation of ice accretion. Section 3 reviews the different CFD modelling and simulation of ice accretion studies reported in the recent literature, including 2D, quasi-3D and 3D simulations. Finally, Section 4 summarizes the principal results and presents the main conclusions.

Computational Models for Ice Accretion
The ice shape and mass prediction use successive airflow calculations, supercooled water droplets trajectories to determine the captured water mass, and heat and mass balance on the blade surface. The latter is also beneficial for the design of icing protection systems [7]. However, icing is not a steady process. As ice builds up, the airfoil geometry changes progressively, which affects the subsequent flow analysis. Therefore, we adopt Unlike the previous icing studies, this paper, in addition to highlighting the problem of wind turbine icing in northern climates, focuses on the actual state-of-the-art of ice modelling and simulation. It provides guidelines to the knowledge needed for modelling the phenomenon and examines the numerical methods used to analyze ice. It provides the best practices and recommendations for the simulation of ice accretion on wind turbines. After reviewing hundreds of references available in the literature on wind turbine icing, investigations have been carried out on the art of modelling ice accretion, from the applicable theory and various mathematical and physical models to the methodological approaches adopted to conduct reliable modelling in this domain. The objective is to present and recommend trustworthy modelling tools that are capable, via numerical simulations, of predicting and quantifying icing and eventually adapting the wind turbines to handle severe winter conditions.
The present paper is structured as follows: Section 2 summarizes the computational approaches employed in the simulation of ice accretion. Section 3 reviews the different CFD modelling and simulation of ice accretion studies reported in the recent literature, including 2D, quasi-3D and 3D simulations. Finally, Section 4 summarizes the principal results and presents the main conclusions.

Computational Models for Ice Accretion
The ice shape and mass prediction use successive airflow calculations, supercooled water droplets trajectories to determine the captured water mass, and heat and mass balance on the blade surface. The latter is also beneficial for the design of icing protection systems [7]. However, icing is not a steady process. As ice builds up, the airfoil geometry changes progressively, which affects the subsequent flow analysis. Therefore, we adopt an approach that accounts for geometry deformation in multiple time steps for more detailed ice form estimation. This requires a finite loop of calculation until we cover the entire estimated period of ice accretion [7,17,18].
Typically, the computational procedure for ice accretion is a time-step iterative process. It requires considering four consecutive modules of calculation: flow field aerodynamic calculations, particle trajectory calculations, thermodynamic analysis and ice geometry generation. All the modelling techniques mentioned above uses partial differential equations that require numerical solutions. This section outlines the mathematical models that describe the modules to be adopted when modelling ice accretion on wind turbine blades' airfoils.

Aerodynamic Modelling
It is essential to simulate the aerodynamic behaviour of clean and iced airfoils to estimate the ice accretion on a wind turbine. Airflow calculations generate the vector field of the continuous fluid phase flow around the airfoil or the wind turbine's blade. The Navier-Stokes equations govern fluid motions. For real fluid flows, where the effect of viscosity is present, these equations are nonlinear. There is no analytical solution for these second-order differential equations. Therefore, numerical solutions are indispensable. For efficient multiphase flow simulations, the most advanced 3D models such as those in the Fluent software use a combination of finite element and finite volume methods to solve the Reynolds-averaged Navier-Stokes (RANS). The last approach makes it possible to analyze real (viscous and compressible) and turbulent flows. This is efficiently used for two and three-dimensional flows [19]. The general form of the momentum equation for an elementary volume for a Newtonian fluid (Navier-Stokes equation) is expressed as: The numerical solution of the RANS equations requires significant time and computational resources. Therefore, assumptions are adopted to simplify these equations, such as considering inviscid fluid, incompressible, irrotational or Newtonian flows. Three simplifying methods have been used in aeronautics to compute the flow field: conformal mapping method, panel method and Eulerian approach [20,21]. The two latter methods are also implemented for wind turbines [22][23][24][25]. Limitations should be noted when applying these simplifying methods [20]. The aerodynamic codes that use the Panel method or Euler's equations (non-viscous fluid) for modelling two-dimensional non-viscous potential flows have obvious limitations over the quality of the results [14].
Alternatively, a combination of the two methods can be used; the panel method is used to facilitate simulation, and the Navier-Stokes equations validate the calculations [12]. Homola, et al. [26] used this technique to estimate the effect of temperature and droplet size on ice accretion on the 5 MW virtual NREL wind turbine blades. The numerical simulations were carried out using a combination of a panel approach and computational fluid dynamics-based techniques. The numerical modelling was conducted in two steps: first, TURBICE (a panel method-based code developed at the Technical Research Center of Finland [25]) was used to reproduce the iced profiles, and then CFD simulations using FLUENT were used to compare the aerodynamic characteristics of the clean and iced airfoils. The iced airfoils generated by TURBICE were used in FLUENT to calculate the iced-airfoil aerodynamic coefficients [26]. Figure 3 illustrates the framework underlying the aerodynamic numerical analysis using the Computational Fluid Dynamics (CFD).

Turbulence Modelling
Turbulence can be expressed as a stochastic process. Flow is turbulent in almost all cases of wind turbine applications, extremely in the presence of ice. Due to the complexity of turbulent flows, there is no analytical theory to predict velocity fluctuations in these flows adequately. It is, therefore, important that turbulence modelling is in concordance with the numerical solution of the flow nonlinear partial differential equations. Depending on statistical tools rather than purely physical methods, turbulence may be expressed in time-averaged equations of fluid flow motion. The Reynolds-Average statistical approach is employed in these cases. For that, CFD is usually expressed as the solution of Reynolds-averaged Navier-Stokes (RANS) equations for turbulence modelling. Turbulence models are discussed in detail in the seventh chapter of Schaffarczyk's book within CFD application for wind turbine aerodynamics [27].
Several models are usually used in CFD simulations to account for turbulence evolution inflows. However, three of them, namely, Spalart-Allmaras, k-ε, and k-ω SST have been widely used in flow simulations. These models, implemented in ANSYS Fluent, refer to the Reynolds-averaged Navier-Stokes (RANS) family of turbulence models [28], have been extensively examined in the literature for both clean and iced wind turbines [29,30].

Turbulence Modelling
Turbulence can be expressed as a stochastic process. Flow is turbulent in almost all cases of wind turbine applications, extremely in the presence of ice. Due to the complexity of turbulent flows, there is no analytical theory to predict velocity fluctuations in these flows adequately. It is, therefore, important that turbulence modelling is in concordance with the numerical solution of the flow nonlinear partial differential equations. Depending on statistical tools rather than purely physical methods, turbulence may be expressed in time-averaged equations of fluid flow motion. The Reynolds-Average statistical approach is employed in these cases. For that, CFD is usually expressed as the solution of Reynoldsaveraged Navier-Stokes (RANS) equations for turbulence modelling. Turbulence models are discussed in detail in the seventh chapter of Schaffarczyk's book within CFD application for wind turbine aerodynamics [27].
Several models are usually used in CFD simulations to account for turbulence evolution inflows. However, three of them, namely, Spalart-Allmaras, k-ε, and k-ω SST have been widely used in flow simulations. These models, implemented in ANSYS Fluent, refer to the Reynolds-averaged Navier-Stokes (RANS) family of turbulence models [28], have been extensively examined in the literature for both clean and iced wind turbines [29,30].

SPALART-ALLMARAS Turbulence Model
The Spalart-Allmaras is a one-equation model. This additional equation is used to model the turbulence viscosity transport. Initially developed for aerospace applications, this turbulence model is widely used for wind turbines due to its compromise between acceptable computational cost and the required accuracy for simulating the turbulent flow [31,32].
When simulating ice accretions on conductors, wind turbines and aircraft, some authors [32,33] pointed out that the Spalart-Allmaras model is the best performing. Based on icing simulation studies carried out for aeronautical purposes using both the Spalart-Allmaras and k-ω SST turbulence models, Mortensen [34] found that the solution using the Spalart-Allmaras model converges faster and more easily than that with the k-ω SST model. Makkonen, et al. [33] stated that Spalart-Allmaras as a one-equation turbulence model is easy to solve and is suitable for simulating the airflow during ice accretion on conductors, wind turbines and aircraft.

k-ε Turbulence Model
The k-ε is a two-equation turbulence model. The two equations are based on the turbulent kinetic energy (k) and the dissipation rate (ε). The simplest, the standard k-ε model, proposed by Launder and Spalding [30], assumes that the flow is fully turbulent and, hence, gives better results for ice accretion simulation. As mentioned previously, when ice is accreted on wind turbine blades, the aerodynamic form of airfoils deforms drastically. This leads to more turbulence in the boundary layer on the suction side of the airfoil. A more developed k-ε turbulence model, RNG k-ε, based on renormalization group theory [35], has correction terms for swirling flow, low Reynolds number flow, and flow with high-velocity gradients. Therefore, the k-ε turbulence model is used as a compromise between acceptable computational cost and the required accuracy in simulating the fully turbulent flow [26]. Villalpando, et al. [29] used the RNG k-ε in a study intended for estimating the aerodynamic characteristics of wind turbines blade sections operating under icing conditions. They found, compared to the two other tested models (k-ω SST and Spalart-Allmaras), that the RNG k-ε model gives the highest prediction of lift at maximal lift angle, which is also affected by ice accretion. The model also predicts the smallest recirculation zone and the separation point that is closest to the trailing edge. Villalpando, et al. [29] stated that this model was adequate for simulating icing on wind turbines.

k-ω SST (Shear Stress Transport) Turbulence Model
This model, originally introduced by F. R. Menter [36], combines the standard k-ω model and the transformed k-ε model. It combines the specific advantage of each of these two models in different regions of flow; the standard k-ω model is activated near the blade's surface and the k-ε model away from the wall [30]. In modelling a wind turbine blade, the k-ω SST offers advantages related to studying the laminar-turbulent transition of the boundary layer at high angles of attack [37]. Several studies have used it to account for turbulence in wind turbine aerodynamic calculations [27].
Regarding the effect of icing on the aerodynamic simulation of wind turbines, the k-ω SST model has been examined in several studies for various types of iced airfoils. Villalpando [38] used the k-ω SST model to study flow around an iced wind turbine NACA 63-415 airfoil using the commercial package FLUENT [38]. The k-ω SST model was found to be more accurate to account for the boundary limit separation of the iced airfoils, especially for simulating the operation under high angles of attack [39]. It also can handle the recirculation zones and accurately predict flow separation [29]. The k-ω SST model has been shown to work best for flows with strong adverse pressure gradients [40,41], being able to describe the generation of specific vortices at the trailing and leading edges [42]. The effect of turbulence models on the aerodynamic performance of wind turbines has been discussed in several studies [29,30,37,43]. Pedersen and Yin [44] conducted an icing simulation study on a wind turbine airfoil using the three mentioned turbulence models. Lift and drag coefficients were the variables of interest. The comparisons were made with experimental data and information available from the literature. They concluded that the k-ω SST is the most appropriate turbulence model for icing simulations for wind turbines. Some authors carried out their studies by implementing two turbulence models. For example, Hudecz, et al. [32] used the Sparlart-Allmaras model during the first 500 iterations, providing an initial guess for the k-ω SST turbulence model to achieve more stable convergence.

Modelling of Surface Roughness
Roughness is a critical aspect of icing software accuracy. Blades' surfaces have a high sensitivity to roughness since it affects the boundary layer transition resulting in flow separation [23,45]. Any formation of ice on the airfoil can affect the turbine's energy output [16], and subsequently, the energy required by the icing protection systems. It is even worthwhile to account for roughness effects, which vary throughout the evolution of ice shape growth, for every temporary stage of ice accretion estimation [23].
Several aircraft icing studies concluded that the roughness height is an important parameter to be integrated into wind turbine icing modelling since it affects the convective heat transfer coefficient estimation. Even with a thin film of ice, surface roughness can considerably increase the local heat transfer [46], which is the critical parameter of heat balance analysis [16,47]. The effect of surface roughness on the heat transfer coefficient has been studied theoretically and numerically by Makkonen [48] and Szilder and Lozowski [49].
Among the most common correlations used to estimate the icing surface roughness over wind turbine blades is the Shin, et al. [50] sand-grain roughness model initially developed for aeronautics. The empirical correlation depends on the Shin and Bond formula, which estimates the small-scale surface roughness height k/c (mm) as a function of liquid water content (LWC) (g/m 3 ), median volume diameter (MVD) (µm), airfoil chord length(c) (m), static temperature(T) ( • C) and the relative wind speed(V) (m/s) considered in the simulation. Based on the validation data available in the literature, a sand-grain roughness value of 5 mm has been determined to be an appropriate setting for most icing calculations [51]. Furthermore, the advanced commercial icing software runs implicit surface sand-grain roughness distribution calculations using the parameters considered for simulation. For example, FENSAP-ICE can integrate constant and variable sand-grain roughness distribution in its calculation [51].

Multiphase Modelling of Droplet Trajectories
The objective of the computation of supercooled water droplet trajectories is to figure out the ability of an airfoil to capture water droplets present in the flow. It is demonstrated through the graphical representation of the local collection efficiency along the airfoil curvilinear (see Figure 4). The local collection efficiency (β) represents the ratio of massflux of the impinging droplets to the mass-flux in the free stream. The collection efficiency is determined based on the velocity vector, and it allows for the estimation of the mass flow rate captured by the surface [52,53].  For 2D problem calculations, as in the case of an airfoil, two impingement limits have to be determined for particle trajectory. For three-dimensional flow, the impingement limits may vary spanwise along the surface of the blade. When considering the rotation of the airfoil, the problem is no longer 2D, and it requires the adoption of four trajectories rather than two [7,54]. Some numerical codes consider the 3-D approach for the water impingement process (LEWICE, ONERA, FENSAP-ICE). It should be noted that droplets diameter affects the impingement limits determination. Median volume diameter (MVD) distribution can be used to calculate collision efficiency, showing a good approximation for a droplet spectrum [26]. When MVD approximation is used to calculate local collision efficiency, the impingement limits have to be established for each droplet size, and the impingement limits of the largest droplet diameter in the distribution are considered [7]. Figure 5 demonstrates an overview of droplet trajectory module calculations.
Ice accretion is a multiphase flow problem. It involves the combination of air as a continuous phase and supercooled water droplets as a dispersed phase. Water droplet trajectory calculations make it possible to obtain the local collection efficiency for complex geometric shapes, using two common approaches: the Eulerian approach and the Lagrangian approach [12]. The two approaches describe fluid motion in multiphase flows that include particles or droplets in two different ways [55]. However, the local collection efficiency can be obtained with satisfactory results using either Lagrangian or Eulerian approaches. The advantage of each of the two approaches is linked to their precision. Otherwise, the efficiency should therefore be verified since the Lagrangian approach comprises a relatively high operating cost but possibly of the same order as the Eulerian approach for the same precision (compared to the mass flow hitting the surface). For 2D problem calculations, as in the case of an airfoil, two impingement limits have to be determined for particle trajectory. For three-dimensional flow, the impingement limits may vary spanwise along the surface of the blade. When considering the rotation of the airfoil, the problem is no longer 2D, and it requires the adoption of four trajectories rather than two [7,55]. Some numerical codes consider the 3-D approach for the water impingement process (LEWICE, ONERA, FENSAP-ICE). It should be noted that droplets diameter affects the impingement limits determination. Median volume diameter (MVD) distribution can be used to calculate collision efficiency, showing a good approximation for a droplet spectrum [26]. When MVD approximation is used to calculate local collision efficiency, the impingement limits have to be established for each droplet size, and the impingement limits of the largest droplet diameter in the distribution are considered [7]. Figure 5 demonstrates an overview of droplet trajectory module calculations.
Ice accretion is a multiphase flow problem. It involves the combination of air as a continuous phase and supercooled water droplets as a dispersed phase. Water droplet trajectory calculations make it possible to obtain the local collection efficiency for complex geometric shapes, using two common approaches: the Eulerian approach and the Lagrangian approach [12]. The two approaches describe fluid motion in multiphase flows that include particles or droplets in two different ways [56]. However, the local collection efficiency can be obtained with satisfactory results using either Lagrangian or Eulerian approaches. The advantage of each of the two approaches is linked to their precision. Otherwise, the efficiency should therefore be verified since the Lagrangian approach comprises a relatively high operating cost but possibly of the same order as the Eulerian approach for the same precision (compared to the mass flow hitting the surface).

Lagrangian Approach
With this approach, we track the supercooled water droplet trajectories as they cross the mesh. The precision of this method depends on the number of trajectories considered and the temporal discretization of the solution diagram. The velocity field is also an important factor in the accuracy of this method, which is directly related to the choice of the flow solver. Therefore, the Lagrangian method applies flawlessly to a flow solver with velocity potential, although it could also be applied to a Navier-Stokes or Euler solver. To be precise, the Lagrangian method must have recourse to several trajectories very close to each other, or even better, to call upon interpolation methods. The latter is relatively truthful insofar as the geometry does not offer significant discontinuities and where the flow is therefore locally homogeneous [56].

Eulerian Approach
This approach considers droplets in the air as a continuous phase and uses the concept of droplet volume fraction to represent the amount of water within a given control

Lagrangian Approach
With this approach, we track the supercooled water droplet trajectories as they cross the mesh. The precision of this method depends on the number of trajectories considered and the temporal discretization of the solution diagram. The velocity field is also an important factor in the accuracy of this method, which is directly related to the choice of the flow solver. Therefore, the Lagrangian method applies flawlessly to a flow solver with velocity potential, although it could also be applied to a Navier-Stokes or Euler solver. To be precise, the Lagrangian method must have recourse to several trajectories very close to each other, or even better, to call upon interpolation methods. The latter is relatively truthful insofar as the geometry does not offer significant discontinuities and where the flow is therefore locally homogeneous [57].

Eulerian Approach
This approach considers droplets in the air as a continuous phase and uses the concept of droplet volume fraction to represent the amount of water within a given control volume. The droplets velocity and the volume fraction of water are only computed at the nodes using the same grid and numerical techniques [57].
The Eulerian approach consists of establishing mass balances on the adjacent control volumes by solving mass, momentum, and energy conservation equations and resolving the Navier-Stokes equation for the fluid. The Eulerian approach needs a mesh in space and a solution for the velocity field on this mesh. The precision depends on the mesh quality and refinement, which are often dictated by the flow problem and not by the droplet trajectories or ice shapes. In addition, it must solve transport equations for each control volume over the entire domain, which adds to the simulation cost. On the other hand, there is no need to perform scanning to calculate trajectories or monitor whether there is an impact [7,58,59]. A clear advantage of the Eulerian approach is its evaluation of the mass flow, which is only approximated by the spatial discretization.
The Eulerian approach is more appropriate for determining the mass of water droplets that impinge on the surface [7]. This is because it defines a clear idea of the entire flow. Within the Eulerian approach, the local collection efficiency is calculated by the inner product of the normal surface vector and the velocity vector of the droplets times the droplet volume fraction [60]. Pedersen and Yin [44] also concluded that the Eulerian multiphase model is an appropriate sub-model for multiphase flow icing calculations for wind turbines [44].

Thermodynamic Modelling of Ice Accretion
This model is essential for analyzing the thermodynamic characteristics of the icing process on wind turbine blades during the impingement of the water droplets by considering the mass and energy balance on the blade surface. As previously stated in the theoretical concepts, the supercooled water droplets may not automatically freeze on impact like with rime ice accretion. Nevertheless, it may leave a fraction of the liquid water running back along the blade's surface to freeze downstream due to the insufficient heat transfer. This process is strongly related to mass and heat transfer, surface roughness, skin friction, and other relevant domains.
Based on the mass flux and heat balance, the freezing fraction of the incoming water droplets for a control volume can be calculated. Along with the droplet impingement computation, the amount and growth of ice formed in that control volume can be determined. The ambient temperature is critical for determining the type of surface involved: dry, wet, or liquid, and therefore the energy balance. For a temperature near the freezing point, the blade surface is fairly liquid, whereas, for glaze ice conditions, the surface is almost entirely wet. For rime ice conditions, the blade surface is mostly or entirely dry [12].
The thermodynamic analysis is traditionally based on the first-order differential equations of mass conservation and heat transfer developed by Messinger [61] and on the work of Lozowski [61][62][63][64]. The model of Messinger was initially developed for the measurements of LWC and to calculate the energy required for icing protection. In addition, Messinger introduced the concept of freezing fraction, as the fraction of impinging water that freezes within the impingement region. However, the model gives moderately good results for rime ice while the results are not enough satisfying for the glaze [7]. Figure 6 illustrates the thermodynamic module calculations based on the Messinger model for mass and energy balances.
For most ice simulation models, the liquid phase is not adequately represented, which is a simplification that limits the prediction quality of these models. However, this phase is always present, even in the dry mode, for very short moments, and it dominates the form of ice accretion [12]. The significance of the wet state phenomena is generally to assume that all unfrozen liquid, generated on a surface element, is completely entrained during a time increment and flows to the next surface element. When the supercooled water droplets do not freeze instantly upon hitting the object, the water droplets caught by the object's surface will coalesce and form surface drops. These drops, which are larger in diameter than the water droplets, make the surface element rough. It should also be noted that, without detailed simulation of the liquid phase, the local roughness, ice density, residual liquid water, chipping and splashing of supercooled water droplets are generally estimated using empirical correlations [12].
As an alternative for CFD, empirical models have been developed to quickly simulate the ice accretion on structures under different meteorological conditions. For example, Makkonen model, referenced in the ISO 12494 standard, is commonly used for ice accretion on cylinders [10,65]. This model, based on three fractions: Collision efficiency, sticking efficiency and accretion efficiency [66,67], was later updated to include a more detailed treatment of wet snow growth to have more flexibility for a wide range of icing problems [68]. noted that, without detailed simulation of the liquid phase, the local roughness, ice density, residual liquid water, chipping and splashing of supercooled water droplets are generally estimated using empirical correlations [12]. As an alternative for CFD, empirical models have been developed to quickly simulate the ice accretion on structures under different meteorological conditions. For example, Makkonen model, referenced in the ISO 12494 standard, is commonly used for ice accretion on cylinders [10,64]. This model, based on three fractions: Collision efficiency, sticking efficiency and accretion efficiency [65,66], was later updated to include a more detailed treatment of wet snow growth to have more flexibility for a wide range of icing problems [67].

Research Survey on Modelling and Simulation of Ice Accretion on Wind Turbines
Comparative information of the different studies of ice accretion on wind turbines is provided in a tabular format. Table 1 organizes information chronologically and provides the roughness model, the turbulence model, the principal parameters calculated (output

Research Survey on Modelling and Simulation of Ice Accretion on Wind Turbines
Comparative information of the different studies of ice accretion on wind turbines is provided in a tabular format. Table 1 organizes information chronologically and provides the roughness model, the turbulence model, the principal parameters calculated (output parameters), and if the authors carried out a process of validation (experimental and numerical). This section analyzes the different modelling approaches by evaluating their suitability and adequacy for wind turbine icing simulation to identify the most adapted techniques and methodologies.

Research on Surface Roughness Modelling
The roughness effect on WT is discussed thoroughly in the Ph.D. thesis of Sagol [16]. However, many other authors have focused their studies on the effect that roughness has on other parameters. For example, Hansman and Turnock [105] conducted a survey of the behaviour of the fluid on a rough cylinder to observe the laminar-turbulent transition of the flow and also to evaluate the heat transfer coefficient. Hansman Jr and Turnock [106] demonstrated that surface tension is probably the main factor responsible for the formation of droplets on wing surfaces. Broeren and Bragg [107] found that ice roughness is generally more prominent than the local boundary layer. Increasing height increases the aerodynamic effect, having a significant impact on separation downstream. Achenbach [108] found a strong dependence between the local convective heat transfer, the ice roughness, and the boundary layer behaviour. General roughness effects are also discussed in Timmer [109]. Etemaddar, et al. [81] mentioned a study where they found that icing mostly affects the roughness on the first 25% of the airfoil chord. However, roughness does not have determinant factors for analysis over the ice-affected surfaces. Therefore, roughness scale parameters estimation is almost empirical [46].
It is difficult to well estimate the roughness as the film of ice is evolving over time. Most of the empirical models have been developed for aeronautics. The methodology consists of running many icing tunnel experiments to create iced airfoil shapes, and then the ice shapes are reproduced through CFD to came up with correlations. Most of these correlations are nondimensionalized for the thickness-to-chord ratio of aircraft airfoils. They are not tested for accretion over other shapes of the structure. It is therefore essential to predict the roughness distribution over the blade for which icing is being simulated. According to Fortin, et al. [62], empirical roughness models developed for aeronautics are not well adapted to wind turbine blades due to differences in operational and meteorological conditions. The work of Fortin [12] led to the development of a surface roughness model suitable for the wet type of ice that can be adapted to wind turbines icing [62,110].
For wet growth icing, beads grow on the surface affecting the roughness, and the empirical models become unrealistic. In this case, FENSAP-ICE is equipped with an analytical roughness prediction code that can show the roughness distribution over any structure in space and in time. In this software, the use of the "beading model" (which integrates roughness change throughout the evolution of ice shape) together with the "multi-shot approach" is beneficial to account for roughness distribution variation, considering moving and freezing beads, throughout the evolution of ice shape. The initial value is calculated depending on one of the roughness models included in FENSAP-ICE [51]. In a multi-shot process, roughness data can be transferred to the flow solver to be considered for calculating shear stresses and heat fluxes. After each shot, the grid is updated to conform to the new ice shape, providing the surface roughness of ice formation if the beading model is enabled [51]. To examine the effect of considering roughness throughout the phases of growing ice simulation, we conducted a Multi-shot simulation of glaze ice on a NACA 64-618 airfoil using the Beading model in FENASP-ICE (integrating roughness change throughout the evolution of ice shape). The simulation was carried out with three shots of 20 min each for a total of 60 min of ice accretion at T = −9 • C. The results on the iced form generated for every shot are illustrated in Figure 7, showing more accuracy and details on the boundary of the formed ice. As mentioned previously, roughness increases the convective heat transfer and the local collection efficiency, which increases the growth of ice accretion in this region [61]. The increase in the collection efficiency due to beading roughness has never been considered in the calculations of local collection efficiency using CFD. However, corrections for local collection efficiency are proposed by Yoon and Yee [59] based on the roughness element concept.

Research on Droplets Trajectory Modelling
Traditionally, most ice accretion codes have used the Lagrangian approach, while the recently developed ones use the Eulerian approach. As an example, LEWICE and ANSYS-FENSAP-ICE, are representative codes based on the Lagrangian and Eulerian approaches, respectively. Homola, et al. [31] simulated the two-phase flow using the Eulerian-Eulerian model in FENSAP-ICE. The main advantage of using this approach is that the same mesh can be used for multiphase flow calculations and ice geometry (Homola,et al. [31], Virk, et al. [75]). Pedersen and Yin [44] used the Eulerian multiphase model droplet modelling as it is "easier" to estimate collection efficiency in a 3D-based simulation. The STAR-CCM+ code from Siemens Simcenter uses a combination of hybrid multiphase models such as Eulerian Multiphase (EMP), Dispersed Multiphase (DMP) and Lagrangian Multiphase (LMP). The latter two were used in recent studies on wind turbine icing, demonstrating accuracy in evaluating the aerodynamic performance drop and the power loss due to icing [39]. Another method for calculating the local collection efficiency of wet ice accretion was proposed by Szlider and Lozowski [48] using morphogenetic models. It is simulated using empirical, boundary-layer or full Navier-Stokes models for local heat transfer [45].

Research on modelling the Electro-Thermal Icing Protection Systems
For wind energy studies, icing prevention is as essential for ice modelling as the ice As mentioned previously, roughness increases the convective heat transfer and the local collection efficiency, which increases the growth of ice accretion in this region [62]. The increase in the collection efficiency due to beading roughness has never been considered in the calculations of local collection efficiency using CFD. However, corrections for local collection efficiency are proposed by Yoon and Yee [60] based on the roughness element concept.

Research on Droplets Trajectory Modelling
Traditionally, most ice accretion codes have used the Lagrangian approach, while the recently developed ones use the Eulerian approach. As an example, LEWICE and ANSYS-FENSAP-ICE, are representative codes based on the Lagrangian and Eulerian approaches, respectively. Homola, et al. [31] simulated the two-phase flow using the Eulerian-Eulerian model in FENSAP-ICE. The main advantage of using this approach is that the same mesh can be used for multiphase flow calculations and ice geometry (Homola,et al. [31], Virk, et al. [76]). Pedersen and Yin [44] used the Eulerian multiphase model droplet modelling as it is "easier" to estimate collection efficiency in a 3D-based simulation. The STAR-CCM+ code from Siemens Simcenter uses a combination of hybrid multiphase models such as Eulerian Multiphase (EMP), Dispersed Multiphase (DMP) and Lagrangian Multiphase (LMP). The latter two were used in recent studies on wind turbine icing, demonstrating accuracy in evaluating the aerodynamic performance drop and the power loss due to icing [39]. Another method for calculating the local collection efficiency of wet ice accretion was proposed by Szlider and Lozowski [49] using morphogenetic models. It is simulated using empirical, boundary-layer or full Navier-Stokes models for local heat transfer [46].

Research on Modelling the Electro-Thermal Icing Protection Systems
For wind energy studies, icing prevention is as essential for ice modelling as the ice growth on wind turbine blades. A wind turbine icing model is supposed to simulate both ice accretion and ice prevention (active blade heating using anti or de-icing strategies) to provide information for investment analysis and optimize icing solutions [111]. Optimizing wind turbines' operation in icing conditions essentially depends on optimizing the operation of the Icing Protection Systems (IPS), given that these systems are energy-consuming and still need more development to fulfill the desired maturity for wind turbines [112]. Modelling is an essential tool to estimate the heat and the operation time required from the IPS to obtain a fully ice-cleaned surface of the blade for a given icing event during a targeted icing duration. The design and optimization of the electro-thermal IPS imply conducting various scenarios of de-icing or anti-icing over the different surface areas of the blade. This should define the protected zones, heat load and critical points to optimize the process by lowering costs and energy consumption.
Calculation of ice accretion depends primarily on thermodynamics, which determines the amount of energy required in the form of heat for a system to pass from one state to another, including phase change. In contrast, heat transfer considers the mechanisms and the rate at which heat is exchanged [113]. For IPS modelling, a module for heat transfer is conjugated in the calculations. This includes the calculations of the heat transfer coefficient using a boundary layer calculation considering the roughness of the wall due to the presence of ice. Some commercial codes are equipped with this additional module to conduct simulations on the hot air or electro-thermal icing protection systems to estimate the energy consumption required for anti-or de-icing once the ice accretion is estimated. For example, LEWICE integrates a thermal anti-icing model for the calculation of the thermal energy required to prevent ice formation, available for two anti-icing modes of prevention: evaporative and running-wet conditions [83]. FENSAP-ICE also has a particular module (CHT3D) that incorporates heat flux distributions in icing calculations to estimate the required total heat load for the hot air or electro-thermal IPS for the two modes of prevention: evaporative and running-wet conditions [51].
Performance degradation due to Icing and heat flux simulation for different anti-icing scenarios are presented in a simulation study on wind turbine blades using FENSAP-ICE [79,80]. This showed the total power required for every considered surface coverage and power distribution for four icing cases. It is worth noting that for initial de-icing computations, the surface roughness should be considered in calculations as it is already determined during ice accretion. However, for anti-icing simulation, uncontaminated surfaces can be considered (with no roughness) since they are expected to remain free of ice [51].
Modelling of the electro-thermal Icing Protection Systems for wind turbines presents many challenges. They are presented and discussed, along with the mathematical models used in IPS modelling, in a published study by Fakorede, 2018 [114] on the modelling of heat and mass transfer during the anti/de-icing process for wind turbine blades. The thermodynamic models used to calculate de/anti-icing systems are also presented in a study by Fortin and Perron [18] on an airfoil NACA 63-415. The paper also includes description of atmospheric icing, icing parameters and the different types of ice accretion, in addition to a review of the thermodynamic models used in ice accretion calculations for wind turbines and the numerical methods used for water droplet trajectory calculations and aerodynamic performance degradation [18]. The various published studies cited on wind turbine icing simulation were mainly up to quite recently on ice growth simulation. Of those, few studies integrate the modelling of electro-thermal IPS in their calculations [111]. In the last decade, ice accretion research has intensified on several techniques for ice prevention solutions. Comparisons between mitigation techniques, their executive design, verification and certification, advantages and disadvantages were systematically discussed in several studies [7,11,15,[115][116][117][118] and Chapter 5, pp. 251-324 in [7]. A significant report by the Technical Research Centre of Finland (VTT) published in 2018 analyzed the performance and the maturity of wind turbines equipped with Ice Protection Systems for various icing climates. The report discussed IPS maturity from the viewpoint of availability and heating durations deviation analysis [112].
Numerical attempts have been made to make virtual validations and test multiple combinations and strategies before conducting costly experimental validation. The numerical models are mainly based on computational fluids dynamics, heat transfer and optimization theories. References [8,18,25,114,[119][120][121][122][123][124][125][126] are examples of these numerical studies on thermal de-icing and anti-icing systems. The recent tendency of wind energy IPS research is related to the pulse electro-thermal de-icing technique [114]. Additional parameters should be considered in an IPS study, such as thermal flow, operational mode, system control, icing detection, the effectiveness of the protection, operating and investment cost, and risk of ice ejection [114].

Conclusions
In this paper, a comprehensive review of wind turbine icing modelling is presented. The commonly applied approaches and the recent advances have been reported based on a systematic and concise literature review. The paper presents a panorama of the current modelling and simulation studies dealing with different aspects of ice accretion on wind turbines. The review of numerous CFD-based published studies on wind turbine icing during the last decade leads to the following considerations: The most common approach uses Computational Fluid Dynamics (CFD). It treats multi-physics aspects including aerothermodynamics, thermodynamics, multiphase flows, heat and mass transfer and optimization theories. Several papers combine the CFD flow and icing simulations with the Blade Element Momentum Theory to generate the power curves for iced-up wind turbines. Some of these attempts were implemented using in-flight icing codes. In contrast, others were special-purpose codes, while most of them were merely methodologies developed using commercial computer-aided engineering software As iced airfoils deform drastically, flow separation is very likely on the suction side. Therefore, turbulence modelling should be carefully investigated for accuracy concerns, especially when closer to the stall angle. Analyzing airfoils beyond stall leads to erroneous results due to flow separation, spanwise flow, and vortices. Therefore, the two-dimensional analysis will no longer be appropriate. For the same reason, software using the panel method is not recommended for modelling iced airfoils.
Another crucial factor to account for in simulations is surface roughness modelling. For the accuracy of ice shape boundaries, roughness should be considered in every interval possible of geometric deformation during the development of ice shape in simulation (combining the multi-shot simulation with the Beading model in FENASP-ICE as demonstrated in Figure 7). The distribution of roughness in time and space over blade surfaces due to ice accretion significantly affects the ice accretion and the aerodynamic characteristics.To date, other than what was cited previously, there is scarce evidence of studies that directly attack the measuring of roughness distribution. The majority of models used the Shin et al. roughness model derived from aeronautics. Automatic meshing and mesh optimization should be emphasized.
As previously mentioned, experimental studies require specialized measuring elements and a special design to analyze the effect of ice on wind turbines. This process is costly and, in some cases, not very accurate. Therefore, little information is available to validate the models proposed in recent years adequately. Thus, the authors generally validated with other numerical studies, other software, experimental data available for specific airfoils (NACA 0012 and S809, mainly), or carried out experimental studies.
It is also noteworthy that almost all the published literature on wind turbine icing studies discussed the icing problem for horizontal-axis wind turbines (HAWTs). On the other hand, very few studies evaluated icing on vertical axis wind turbines (VAWTs) even though VAWT blades are more vulnerable to supercooled water collection due to three-dimensional special impingement. This especially occurs for high rotational speed conditions [127]. However, VAWTs are not common compared to HAWTs, which are more efficient for small-and large-scale power applications.
The control and the optimization studies of the electro-thermal IPS via simulation are poorly presented in the literature, and exclusively in research projects. Most of the research in this domain was in-flight anti/de-icing oriented and adapted to wind turbines. Anti/deicing models should incorporate an ice accretion model for better performance estimations.
Most of the research on wind turbine icing modelling was derived from the aeronautical industry, which has come a long way in this area. Following this progress in aeronautics, intensifying the modelling research on wind turbine icing, for example, by categorizing the wind turbine iced airfoils, could boost our comprehension of this phenomenon, towards practical solutions. Following examples from aeronautics, while considering the particularities of wind turbine icing, was our guideline to make this study a unique reference for researchers in this domain. The aim is to suggest appropriate methods, references, and tools to conduct reliable icing modelling for wind turbines. This study can contribute to advancing the research on this phenomenon, common and frequent in Nordic countries, and to adapting wind turbines and optimizing their operation under icing conditions.