A Review on the Current Status of Icing Physics and Mitigation in Aviation

: Icing on an aircraft is the cause of numerous adverse effects on aerodynamic performance. Although the issue was recognized in the 1920s, the icing problem is still an area of ongoing research due to the complexity of the icing phenomena. This review article aims to summarize current research on aircraft icing in two fundamental topics: icing physics and icing mitigation techniques. The icing physics focuses on ﬁxed wings, rotors, and engines severely impacted by icing. The study of engine icing has recently become focused on ice-crystal icing. Icing mitigation techniques reviewed are based on active, passive, and hybrid methods. The active mitigation techniques include those based on thermal and mechanical methods, which are currently in use on aircraft. The passive mitigation techniques discussed are based on current ongoing studies in chemical coatings. The hybrid mitigation technique is reviewed as a combination of the thermal method (active) and chemical coating (passive) to lower energy consumption.


Introduction
Aircraft icing has always been a hazardous issue since it was recognized in the 1920s because it causes degradation in the aerodynamic performance and malfunction of essential flight instruments [1].The primary focus of aircraft icing research has been on ice accretion on wings for its serious adverse consequences.When icing occurs on a wing, the change in airfoil shape results in a decrease in lift and an increase in drag, leading to potentially fatal accidents [2].Historical icing studies were documented several times in the literature [1][2][3].The early stages of icing research were carried out mainly by experiments geared toward understanding icing physics and the development of icing mitigation techniques [4].In the 1970s, the development of numerical simulations for aircraft icing began growing.The motivation of the simulation was reducing the cost and lowering the risk of accidents during the certification process for icing mitigation devices.Since then, experimental and numerical investigation have been used side by side to establish more accurate icing simulation tools.
From a meteorological perspective, initial icing research focused on icing caused by small water droplets up to 40 µm in diameter, often referred to in Appendix C of Part 25 of Title 14 in the U.S. Code of Federal Regulations [5].However, an ATR-72 accident in Roselawn in 1994 [6] revealed the need for research expertise beyond the water-droplet-size envelope provided in Appendix C. The National Transportation Safety Board concluded that the accident was caused by supercooled water droplets, which were up to 2000 µm in diameter.Therefore, the FAA introduced new regulations as Appendix O, which describes meteorological conditions with potential for icing with water droplet diameter up to approximately 2000 µm in diameter.Since the 1990s, icing caused by water droplets larger than 40 µm in diameter has been referred to as supercooled large droplet (SLD) icing.It has remained a challenging topic for aircraft icing researchers.In addition to SLD icing, the 1990s saw a significant increase in the interest in jet engine icing [7,8].Engine icing is often referred to as ice crystal icing since it is considered to be caused by ice crystals rather than the supercooled droplets that cause icing on wings.Meteorological conditions for the ice crystal icing were documented as Appendix D of Part 33 of Title 14 in the U.S. Code of Federal Regulations [9] and Appendix P of European Union Aviation Safety Agency (EASA) CS-25 [10].For the ice crystal icing, knowledge of the underlying physics is currently limited.
There are two different concepts for icing mitigation: anti-icing and deicing [4].Antiicing techniques are designed to completely prevent ice from forming on an aircraft surface, while deicing techniques remove the accreted ice before it causes significant adverse effects [4].Nowadays, several mitigation methods are available and used for aircraft.However, there is still considerable research towards optimizing the established methods and developing new icing mitigation systems.
It is beneficial to review the forefront of icing physics and icing mitigation techniques to clarify the remaining challenges for icing research.From time to time, aircraft icing research is reviewed by several researchers.Previous reviews discussed the icing physics on fixed-wing aircraft and currently in-use icing mitigation techniques.Gent in 2000 [11] and Cebeci in 2003 [12] focused on the status of calculation techniques for ice accretion prediction, aerodynamic performance degradation, and ice mitigation system performance.Lynch and Khodadoust [13] reviewed test results for aerodynamic performance and control degradation, taking into account several types of ice accretions on a wing.Bragg et al. [3] provided the history of icing research from the beginning and reviewed iced airfoil aerodynamics based on flow field physics.Cao et al. reviewed icing research in general, from meteorological icing conditions to the analysis of actual flight accidents [14].In terms of the icing mitigation techniques, Thomas et al. [4] introduced the available mitigation methods at the time of their publication.
This review article aims to summarize the current knowledge on icing physics and mitigation techniques for aircraft.The research status of icing physics is provided in three parts.The first part focuses on the icing on a fixed wing, which has a long history with numerous mitigation techniques.The second part is given for rotor icing, which has more complicated icing physics than a fixed wing.The last part presents the engine icing problem.The icing mitigation technique begins with discussing active mitigation methods such as pneumatic boots and hot-air anti-icing systems currently used on commercial aircraft.The review then moves toward discussing passive mitigation techniques currently under development.The hybrid mitigation technique combining the active and passive mitigation methods is introduced as a new approach.

Icing Physics
Icing research began with studying icing on fixed wings which can be considered as a two-dimensional phenomenon.On rotors, the three-dimensional effects of the flow field around the airfoil make the icing phenomena more complicated than that of the fixed wing.These two research topics have long histories, and so several simulation codes for ice accretion prediction are available.Additionally, recent experimental research has been conducted to improve the accuracy of the predicted ice accretion shape of those simulation codes.Therefore, both experimental and numerical investigations are presented along with the improvements in the icing simulation techniques.
Engine icing is discussed in the last part of this section.This topic has garnered significant attention since the 2000s.Compared to icing on airfoils, the knowledge of engine icing is limited [15].Recent efforts have been devoted to a better understanding of the icing phenomenon itself.Current findings are summarized to accelerate future exploration.

Fixed Wing Icing
The impact of supercooled water droplets causes the icing on a fixed wing.Depending on meteorological conditions such as temperature, airspeed, and the size of water droplets, ice accretion is formed as rime, glaze, and mixed ice.Rime ice is likely to be created when the temperature is low enough for the entire droplets to freeze entirely upon impact.On the other hand, when water droplets freeze only partially upon the impact (for instance, when the temperature is lower than, yet close to, 0 • C), glaze ice is formed on the surface.Mixed refers to the coexistence of rime and glaze ice.[14] These ice formations can be seen from icing caused by supercooled water droplets up to 40 µm in diameter.However, if the droplet diameter exceeds 40 µm, known as supercooled large droplets (SLDs), further study is necessary to understand the features of ice accretion.In the following, recent icing research is presented after introducing current icing simulation schemes.
Icing simulation tools became available in the 1990s [16].The initial function of those tools was to predict the ice shape on a two-dimensional airfoil from the icing environmental conditions.Numerous discussions in the literature explain how those tools calculate ice shapes [4,13,17].A conventional icing simulation starts from determining the water collection efficiency, which represents how often water droplets impinge on the airfoil.This calculation often refers to the droplet trajectory calculation [11].In the droplet trajectory calculation, there are two ways to consider the size of water droplets [13].One way is to use a droplet diameter distribution, such as Langmuir D distribution [18].Another considers a single droplet diameter, referred to as median volumetric diameter [13].Once the number of impinging droplets is known, the energy and mass balance within a control volume on the airfoil surface is considered, as illustrated in Figure 1.The energy balance model based on thermodynamics for an unheated surface was proposed by Messinger [19].Mass conservation is considered based on the amount of liquid water in a control volume.The mass coming into a control volume consists of (1) impinging water droplets,  By introducing the freezing fraction, f, into the mass balance, the mass of the ice accretion and water flow to the downstream can be expressed as . .
The freezing fraction represents the number of impinging droplets that freeze in a given control volume compared to the total incoming water mass.It is used to define the type of accreted ice depending on its value.When f is 1.0, the accreted ice is classified as rime ice.In contrast, when f is close to 0, the ice is classified as glaze ice.
For energy conservation, contributing energy terms in the control volume are (1) .Q conv : convective heat, (2) .Q imp : kinetic energy of water droplets, (3) .Q latent : latent heat, and (4) .Q sensible : sensible heat.By taking all these terms, a heat balance can be expressed as . .
Q conv and .
Q imp are defined as follows: . . .
where h conv is the convective heat transfer coefficient, T rec is the airflow recovery temperature, T b is the balance temperature of the control volume, and V imp is the velocity of impinging droplets.For .
Q latent , the latent heat of freezing, .
Q f reeze , and the latent heat of evaporation or sublimation, .
Q es , contribute as follows: .
m es (9) where L f and L es are the latent heat of freezing and evaporation or sublimation per kilogram, respectively.
. Q sensible changes with time and temperature.For the ice accretion process in the control volume, it is assumed that the temperature of incoming water increases to the freezing temperature, T f , followed by the ice and unfrozen water reaching the balance temperature, T b .Since the incoming water is due to droplet impingement and water flow from the adjacent segment in the upstream, the terms associated with heating are .Q s, imp = .m imp C p,w T imp − T f (10) .
Q s,r in = .
m r in C p,w T r in − T f (11) where C p,w is the specific heat of the water.When the ice and unfrozen water reach the balance temperature, .
Finally, sensible heat can be obtained from the following equation: .
By substituting Equations ( 2) and (3) into Equation (4), the freezing fraction, f, and the balance temperature, T b , can be derived.When f is derived for every segment on the airfoil, .m ice can be obtained [20].With considerable efforts devoted to developing icing simulations for a fixed wing, many reports [16,21,22] documented that those simulation results showed good agreements with ice shape from experimental data for rime ice accretions.However, there remains the demand for improvement in glaze ice conditions and supercooled large droplet (SLD) icing, as mentioned in previous review papers [3,11,13].Among the parameters given in Equations ( 1) and (4), it is not fully understood which parameters are not adequately modeled.Although under those circumstances, researchers are working to understand the effect of droplet parameters and surface roughness on the accuracy of aircraft icing simulations.
(1) Droplet Impact The droplet trajectory in conventional icing simulations was based on the assumption that water droplets do not deform and coalesce before the impact, nor do they rebound and splash upon impact [23].However, for SLD icing, since the sizes of water droplets are large, the droplets may deform or even break up due to aerodynamic forces.In this case, the assumptions made in the early stage of aircraft icing simulation no longer hold.There were a large number of investigations on droplet deformation and breaking up in the field, such as fuel injection and aerosol atomization.In the aerospace industry, Tan et al. [24] developed a numerical model on droplet deformation and breakup near the leading edge of an airfoil based on their literature review.They reported that the droplets potentially break up where the pressure gradient is severe.Since this numerical study was carried out with limited experimental data, Vargas and Feo experimentally investigated the droplet deformation and breakup further [25].Specific attention was paid to the value of the critical Weber number for droplet breakup conditions in the air.Here, Weber number is defined as follows: where ρ a is the density of the air, V and d are the velocity and diameter of the droplet, and σ is the surface tension between the air and water.As concluded in [24], the critical Weber number to describe the droplet breakup greatly differed depending on the experimental facilities and relative droplet-air velocity.For example, the current version of LEWICE, the icing simulation provided by NASA, added the empirical relationship reported in [26], where the critical Weber number was 13.While LEWICE treats the droplet trajectories in the Lagrangian formulation, another icing simulation, FENSAP-ICE, implemented the Pilch and Erdman model [27] for droplet breakup in the Eulerian formulation [28,29].The critical Weber number in [27] was defined as follows: We critical = 12 1 + 1.077Oh 1.6  (16)   where Oh is the Ohnesorge number defined as Here, µ d is the droplet viscosity, d is the droplet diameter, ρ d is the droplet density, and σ d is the droplet surface tension.
In addition to the droplet breakup, as mentioned in [30], Papadakis et al. observed that there were small droplets that traveled against the airflow near the leading edge of NACA0012 airfoil in the case of SLDs.Those droplets may bounce back or break up upon impact on the airfoil surface, and thus the assumption that all impinged water droplet stays on the surface does not hold.Since then, droplet impingement dynamics have gained significant attention to correct the deficiency in the assumption.The phenomenon of the droplet bouncing back or breaking up upon the impingement has been referred to as the droplet-wall interaction, and explorations to understand the underlying physics are still ongoing.For droplet splashing, five critical parameters need to be determined: the splashing threshold, splashed droplet diameters, splashed droplet velocities, splashing angles, and the mass of sprayed droplets [31].The fundamentals of the splashing behaviors were investigated mainly outside aerospace fields [32][33][34][35], and there exists a large number of models related to the aforementioned splashing parameters.In this review, splashing models which are used in the current icing simulation tools are highlighted.
The volume of fluid (VOF) or momentum of fluid (MOF) methods are known to simulate droplet splashing dynamics with a high degree of accuracy [36][37][38].These methods are capable of describing the time-changing deformation and breakup of a single droplet during splashing.However, their computational costs are considerably high since the mesh size is smaller than splashed droplets.If they were implemented into the aircraft icing simulation, the overall computational cost exceeds current computational capabilities.Therefore, current icing simulation codes utilize semiempirical models.These models are developed to derive the postimpact droplet properties, such as the number of splashed droplets, velocity, and the size of sprayed droplets, directly from the preimpact droplet properties (Figure 2).Since semiempirical models do not consider the complicated physics when water droplets come in contact with the surface, those models can include the effect of droplet splashing in SLD conditions with relatively low computational cost [30].The function of those semiempirical models in droplet impingement is to correct the mass of incoming water droplets into a control volume in Equation ( 1) by considering the mass loss of water due to splashing [30].To calculate the mass loss, those models are required to provide the properties of the postimpact droplets and the threshold for the occurrence of droplet splashing.Most models use nondimensional numbers such as Weber number, Ohnesorge number, and Reynolds number to deal with variations in velocities and diameters of impinging droplets.Unlike the Weber number used to describe the droplet breakup in the air, the Weber number for splashing is defined as follows: where ρ d is the density of the droplet.Ohnesorge number is the same as Equation (17), and Reynolds number is described as follows: For V in Equations ( 18) and ( 19), the normal component of the velocity of the impinging droplet with respect to the surface, V = u 0 cos θ 0 , is used (Figure 2).These models were developed originally in the spraying industry, and established models were introduced to aircraft icing simulations [39].Most investigations were conducted by single droplet impingement experiments [39].In the case of aircraft icing, water droplets may impinge at speeds up to approximately 100 m/s on a dry or wet surface.The surface roughness and the temperatures of the droplet and the surface may also vary in aircraft icing.The accuracy in the prediction of ice shape from a simulation partially depends on selecting an appropriate model to describe various splashing phenomena.
Mundo et al. [40] provided the required description for the threshold and postimpact droplet properties.This model's distinctive feature is that it considers the viscosity of the impinging droplet for the threshold.They defined the threshold for the splashing phenomenon based on Weber and Reynolds numbers of an impinging droplet as follows: where ρ d is the density of the droplet, d 0 is the diameter of the droplet, u n,o is the normal velocity of the droplet with respect to the impinging surface, σ is the surface tension, and µ d is the viscosity of the droplet.The splashing occurs when K M > 57.7.However, this criterion did not include the effect of surface roughness.For the splashing threshold on a wetted surface, Cossali et al. [41] experimentally determined the following relationship based on Weber and Ohnesorge numbers: where δ is the nondimensional film thickness defined as the ratio between the film thickness and droplet diameter.The determining parameter on the left-hand side in Equation ( 18) is often called the Cossali parameter, K C .Later, Trujillo et al.
[42] reviewed existing splashing models, including the Mundo and Cossali models.They provided a combined splashing model that considers the surface roughness and water film thickness on the wetted surface.
From this model, splashing on a dry surface occurs when where R is the normalized surface roughness.They also reported the correlation in the Cossali parameter between dry and wetted surfaces as follows: While the aforementioned splashing models were derived from experimental results, those investigations were performed at lower velocities than the conditions seen in aircraft icing.When those models were implemented into aircraft icing simulations, researchers had to rely on extrapolated data from experimental results.To make those models more suitable for aircraft icing, researchers modified those models by implementing them into icing simulation codes and comparing the calculated collection efficiency with water impingement results conducted in icing wind tunnels.Wright [43] implemented the Mundo model into LEWICE and modified the description for the splashing threshold.He defined the threshold using the velocity normal to the airfoil surface: where LWC is the liquid water content in the freestream.For the diameters, mass, and velocity of the splashed droplets, he proposed the following descriptions using the Mundo parameter defined in Equation (20): where variables d, m, u, and θ denote the diameter, mass, velocity, and impingement angle of the droplet while subscriptions s, o, t, and n represent the secondary (splashed) droplet, original (before impact) droplet, tangential component, and normal component with respect to the impinging surface, respectively (Figure 2).Honsek [44] proposed an Eulerian formulation of a splashing model based on the Trujillo and Lee model.The proposed model was calibrated by comparing the experimental results from [45] with simulated results obtained by implementing the model into DROP3D, a module calculating the droplet impingement characteristics used in FENSAP.Since this model was established in an Eulerian frame, the modification was made only in the mass loss expression, f m , due to the splashing.
While splashing models developed in the early stage relied on single droplet impingement experiments, most of the models were developed based on comparisons between experimental results of the droplet collection efficiency on airfoils and those produced by aircraft icing simulations in which splashing models were implemented.Although several studies in the literature found that those correlations perform well in aircraft icing simulations, the performance of these models in high-speed flows, such as above 30 m/s, is lacking.Detailed information on high-speed droplet impingement may enable the consideration of additional effects such as the interaction between impinging and splashed droplets.
(2) Surface Roughness The calculation of the freezing fraction, f, is another factor that impacts the accuracy of an icing simulation.Traditional simulations took the control volume approach using the Messinger model shown in Section 2.1 to compute the freezing fraction.This model was developed by considering the heat balance between the air, water droplets, and the airfoil with or without ice.It was further extended by Myers [46] by including air compressibility effects.Of the several parameters related to this heat balance calculation, the surface roughness of the accreted ice has been the most critical factor.This is because surface roughness changes flow characteristics such as boundary layer dynamics that result in changes in the local collection efficiency, convective heat transfer, and consequently, the final ice shape.
Olsen and Walker [47] reported a smooth surface around the leading edge for ice roughness in glaze ice conditions followed by a rough region downstream.Shin [48] conducted qualitative ice roughness measurements and boundary-layer measurements to investigate the effects of the iced surface roughness on the boundary layer.It was confirmed that there was a distribution of irregular surface roughness at the beginning of ice accretion, and the roughness eventually reached a constant value.Anderson and Shin [49] established the correlation between roughness characteristics and nondimensional terms such as the freezing fraction, f, and the accumulation parameter.Those correlations were implemented into simulation codes to predict the maximum roughness height.However, Anderson et al. [50] also revealed that those single parameters were not sufficient to describe the ice roughness in glaze ice conditions.
Recent research has aimed to achieve better roughness prediction schemes [51] and model development to compute the heat transfer coefficient [52][53][54].If topics were to better model actual icing physics, experimental and numerical investigations for heat transfer would be needed to account for the actual ice roughness rather than the sand grain roughness used in conventional icing codes.
The ice roughness measurements were traditionally conducted by digitized pencil tracing, ice shape molding, and multiangle photography [55].Lee et al. developed the laser-based 3D scanning method [56] to resolve the roughness better.This measurement technique leads to developing the statistical characterization of ice roughness by adopting the self-organized map approach [57][58][59].More details on those roughness measurement techniques can be found in [59].

Rotor Icing
Rotorcrafts such as helicopters are operated for a long time at altitudes below 6000 m, where supercooled water droplets may exist in clouds [60].The airfoil of rotorcraft is generally smaller in size than that of a fixed-wing aircraft, which means water droplets are more likely to impinge since the collection efficiency increases when the leading-edge radius decreases, as reported in [61].Therefore, helicopters are more susceptible to icing [62].The icing occurs on the forward front of rotors.This leads to unbalanced loads between the blades, which excites certain vibrational modes.One of the distinct differences between icing on rotors and icing on fixed wings is the existence of the centrifugal forces acting on the accreted ice and water film on the rotor surface [63].Simulating this remains a challenge for current three-dimensional icing simulations.Starting with the explanation of the calculation procedure for the three-dimensional icing simulation, recent advancements in research related to the consideration of centrifugal forces are presented in the following section.
The seriousness of the icing problem on rotors was already recognized in 1977 [64].In the early stages, icing experiments were conducted as in-flight tests where test results heavily depended on meteorological conditions.According to a report from the U.S. Army [65], they had 0.  Based on experimental data obtained in icing test facilities, three-dimensional icing simulation codes have been developed, such as LEWICE 3D [69], FENSAP-ICE [70], and ONERA 3D [71].Even though there are differences among those simulation codes, the calculation schemes are similar to conventional two-dimensional icing simulations.Instead of the two-dimensional flow field calculation, most simulation codes employ the three-dimensional flow field calculation.The calculation of the ice accretion is made on a particular blade section in a way similar to the conventional two-dimensional icing simulation.After calculating the ice shape at each radial position, the flow field around the iced rotor is recalculated in three dimensions.This loop continues until the designated iterations are completed [72].
Due to the fruitful knowledge from two-dimensional icing simulations, current three-dimensional icing simulations have icing prediction capabilities.For instance, Narducci et al. [72] investigated the accuracy of the predicted ice shape by LEWICE 3D by comparing the experimental results obtained during the Helicopter Icing Flight Test program carried out by NASA and the U.S. Army [73].However, it does not take into account the influence of the centrifugal forces [63].Accreted ice on a rotor experiences centrifugal forces, unlike fixed wings.That force is more significant close to the tip, and ice can be detached from the rotor when the centrifugal force overcomes ice adhesion force on the surface.This phenomenon is called ice shedding.It was experimentally confirmed by Busch and Bragg [74].
With the establishment of the AERTS Facility, Brouwers et al. [75] developed the ice shedding model using the ice shape predicted by LEWICE.The model was validated with the propeller icing experiments for a hovering rotor.The ice shedding time and location prediction were within 25% of the experimental results (Figure 4).For more detailed ice shedding mechanisms, Zhan et al. [76,77] developed a threedimensional finite element method code to simulate crack propagation in the accreted ice.The ice shape was generated by FENSAP-ICE.The crack propagation was determined based on the stress acting on the ice, assuming that the direction of the crack was perpendicular to the maximum principal stress direction.While their model provided great insight into the details of ice shedding events, it has not been validated against experimental results since there exists little available data.
In addition to the centrifugal force acting on the accreted ice and causing ice shedding, its effect on the water film is also essential to predict the ice shape accreting on rotors more accurately.In the case of icing on a fixed wing, the motion of the water film is in the chordwise direction due to aerodynamic shear stress.However, in icing on rotors, the centrifugal force on the water film causes motion in the spanwise direction.Zhao et al. [78] developed a numerical method to predict the ice accretion on rotors concerning the centrifugal force acting on the water film.Wang and Zhu [79] also developed a numerical simulation of ice accretion on rotors with this centrifugal force consideration.They concluded that the ice thickness at close to the stagnation point increased while it decreased at the frozen limit.Kelly et al. [80] proposed a numerical approach to calculate the performance degradation of rotors.Their method included the water film motion toward the spanwise direction.When compared with experimental results provided in [75], their method improved the accuracy in terms of predicting the ice shape.
Although the accuracy of the ice accretion prediction on rotors has been improved by considering the centrifugal force, the quantitative evaluation of the accuracy has not been investigated since the available experimental results are limited, especially for the water film movement.Researchers have continued to study these problems over the years.However, continuous efforts are required for further improvement in the rotor icing simulation.

Engine Icing
A large number of unexpected power losses were reported from commercial aircraft at altitudes higher than 22,000 ft near thunderstorms [7].In those regions, water was expected to exist in the form of ice particles.Mason et al. [7] hypothesized that those ice particles contributed to ice accretion inside the engine (Figure 5).Here, engine icing is often called ice crystal icing.In their hypothesis, both liquid and ice particles cause ice accretion inside an engine.Liquid water on the surface inside the engine would slow down ice particles, enabling heat to be removed from the surface.When the surface temperature decreases to the freezing point, ice can begin to accrete.In this situation, only the surface temperature is required to reach the freezing point, which means that this accretion might happen even when the local air temperature is above 0 • C.This unique ice accretion occurs on the static component in the low-pressure compressor and eventually leads to rollback due to the blockage.In addition, this accreted ice can cause damage to compressor blades, rotor blade tip rubs, vibration, surge, and flameout as it is shed from the surface [81].While the development of analytical models is required to support the engine certification process as well as next-generation engine design, this kind of icing physics has not been well understood since there exists a limited number of ground-test facilities that can reproduce the ice particle ingestion into the aircraft engine at cruise conditions [8].To better understand ice crystal icing, numerous research efforts have been conducted based on the hypothetical mechanism proposed in [7].The focus of each research effort is on a specific part of the total process of ice crystal icing.Ongoing research can be seen for ( 1) measurement of the ice crystal icing conditions and (2) a numerical approach.(

1) Measurement on Ice Crystal Icing
To better understand the icing physics of ice crystal icing, environmental conditions in the engine core should be first characterized.In addition to the icing parameters discussed in Sections 2.1 and 2.2, such as air temperature, airspeed, mean volumetric diameter of liquid water droplets (MVD), and liquid water content (LWC), the total water content (TWC) should also be considered for ice crystal icing.TWC is the mass of water droplets and ice particles contained in 1 m 3 of the air.It is expressed as the sum of LWC and the ice water content (IWC), and it is used in the form of the melting ratio, LWC/TWC, since liquid water is the key for ice accretion.Besides TWC, the wet-bulb temperature in the engine is considered in ice crystal icing.It is a function of the temperature, pressure, and humidity of the air [15].As previous works revealed [82][83][84], this temperature determines melting and water evaporation and affects ice accretion behavior in the engine.When the wet-bulb temperature exceeds 0 • C, the melting ratio becomes higher, resulting in slush ice formation.The slushy ice may cause physical damage to compressor blades if it is removed from the surface.In contrast, solid ice forms on the surface when the wet-bulb temperature is below the freezing point.If the ice grows large enough to block the airflow, then flameout may occur [7].
During the establishment of the ground-test facilities for studying ice crystal icing, research challenges seemed to arise from the characterization of TWC conditions.Van Zante et al. [85] measured TWC created in an altitude engine research test facility at NASA's Glenn Propulsion System Laboratory using the isokinetic probe (Figure 6) and multi-hot-wire probe.The isokinetic probe measures TWC by evaporating all hydrometeors it captures [86].It was revealed by Davison et al. [87] that the isokinetic probe can be reliable for high TWC measurement.The multi-hot-wire probe is equipped with several strings.It measures TWC based on the power required to keep the wire temperatures constant.While those two instruments showed similar measurements of LWC in the supercooled liquid cloud condition (no ice crystals present), the multi-hot-wire probe recorded lower LWC values in ice crystal cloud conditions at approximately 60% [9].To achieve more accurate measurements of TWC, analysis based on the existing experimental results and the development of other measurement devices is currently ongoing [88].
(2) Numerical Models of Ice Crystal Icing Physics In addition to researchers' efforts to determine essential factors in icing conditions such as TWC, numerical simulations to predict two-dimensional ice accretion due to ice crystal icing have been developed (GlennIce from NASA [89], IGLOO2D from ONERA [90], FENSAP-ICE from ANSYS [91], and ICICLE from University of Oxford [92]).The primary strategy in the development of ice crystal icing simulations is to introduce the unique features of ice crystal icing into conventional icing simulations of a supercooled water droplet.Each simulation is mainly based on the same process as for the icing simulation on the fixed wing discussed in Section 2.1: (1) airflow around the ice accretion site, (2) droplet and ice particle trajectory [93,94], (3) mass and heat balance, and (4) ice accretion.Most of the current focus of numerical research is on the mass and heat balance calculations [95].Many efforts are made in the extension of the Messinger model (Section 2.1) to deal with the mass flow rate of ice particles [96] as well as the erosion of accreted ice caused by ice particle impingement [97].
The most significant difference between the Messinger model and ice crystal icing simulations is the consideration of the ice particles.While there have been several discussions on the development of models for ice crystal icing [89,98], Trontin and Viledieu [96] provided a straightforward schematic of the water mass rate in a control volume on the ice accretion site, as shown in Figure 7.The following equations are derived based on the conservation of mass rate given in Figure 7:  . .
. m imp (34) where η m is the melting ratio of the impinging ice crystals.While the melting ratio, η m , can be measured, the sticking efficiency, ε s , the fraction of impinging mass that adheres to a surface [99], is left unknown.Trontin and Viledieu [96] developed an empirical model based on the experimental results from Currie et al. [100].They measured the sticking efficiency at the stagnation point of a crowned cylinder over a small value range of the melting ratio at Mach 0.25 and 0.4.
Trontin and Viledieu argued that the sticking efficiency resulted from the competition between the ice accretion due to incoming water droplets or ice crystals and the erosion due to the ice crystals.Since the erosion rate is mostly determined by the tangential velocity of ice crystals, they assume that the sticking efficiency in the glaciated condition can be considered as a function of the melting ratio as follows: They derived this function on the following assumptions: 1.

2.
When only supercooled water droplets exist (η m = 1), the sticking efficiency becomes 1 since all impinging droplets stick to the wall.

3.
In the small η m region in Figure 8, the function can be considered linear such that F(η m ) = K # η m where K # is an adjustable number.

4.
The function is smooth and increases with the melting ratio.
Relation between the melting ratio and sticking efficiency (obtained from [96]).
They expressed the function as polynomials.
Based on Figure 8, K # was calibrated as 2.5.Baumert et al. [101] modified the sticking efficiency in the mixed-phase conditions.In this case, the sticking efficiency of the mixed ice crystal depends on the ratio of liquid to total water content at the wall required for the ice crystal to stick, m stick .
m stick is expressed as where K d is an adjustable parameter and f 1 is the total liquid fraction at the wall given as (39) By introducing the same assumptions made in [8], the function F(m stick ) is assumed to be Note that the variable of this function for the glaciated condition is the melting ratio η m = LWC/TWC while the variable for the glaciated condition is the ratio of liquid to total water content at the wall.The adjustable parameters K d and K # were chosen to be 0.3 and 2.5, respectively.
Since researchers' interest in ice crystal icing started growing recently, knowledge on fundamental physics is still limited, both in terms of the full-size engine scale and the ice crystal scale.As of now, it is clear that a detailed description of ice crystal impingement behavior such as breaking, eroding, and sticking on the clear/iced surface may improve the fidelity of a simulation.

Mitigation Techniques
Along with understanding the fundamental icing physics, the development of icing mitigation techniques is vital to improving airworthiness in icing conditions.Since such activities began during World War II, several methods were introduced to mitigate the adverse effects of icing [4].Those mitigation techniques can be classified as either active or passive mitigations.The active mitigation technique uses energy from an external system, while the passive mitigation technique utilizes the physical or chemical properties to combat icing [102].In addition, there are two distinct concepts in icing mitigation techniques: antiicing and deicing.The anti-icing techniques ultimately prevent ice accretion, while the deicing techniques remove ice deposits accreted on the surface [103].

Active Mitigation
(1) Thermal Method The thermal method refers to the mitigation techniques that use heat to increase the surface temperature to prevent icing (anti-icing) or melt the accreted ice (deicing).The heat can be provided from either exhaust gas from the engine (hot air) or an electrical heater embedded underneath the surface of the wing.The thermal method is the most widely used method of icing mitigation for wings [104].While it ensures a safe flight under icing conditions, energy consumption has been a major topic of discussion.For the bleed air system, 2.5~5% of the core engine mass flow is required according to [105,106].This energy extraction contributes to the increase in fuel consumption.Therefore, lower energy consumption is desirable.The electrothermal method is used for fixed wings and propellers and helicopter rotors where the bleed air system cannot be installed.However, since electrical power is provided by onboard generators and thus limited during the flight [107], a reduction in energy consumption is necessary.
To optimize these devices, the required energy to mitigate icing (anti-/deicing) under given icing conditions must be understood.Some approaches to categorizing the parameters affect the energy requirement for icing mitigation [108,109].However, they are similar to the following:

•
Structural parameters of the wing.

•
Types of materials, wing skin thickness, leading-edge geometry, etc.
Arguments regarding the optimization of the thermal mitigation technique have been developed to focus on the surface temperature distribution of the mitigation site.The surface temperature can be investigated experimentally [110] with icing wind tunnel tests or numerically by utilizing the aircraft icing simulation discussed in Section 2.1.Due to the relatively high cost and limitations in experimental conditions, which cannot cover the icing conditions defined in [5], numerical investigations are often selected [104].
Since aircraft icing simulations have been capable of considering the heat provided from the electrothermal mitigation [69], the piccolo tube parameters such as hole spacing and Reynolds number of impinging hot air extracted from the engine are added to deal with the hot-air icing mitigation system in the form of correlations [111,112].Along with the improvement of aircraft icing simulations, a parameter sensitivity analysis has been conducted.Pourbagian and Habashi conducted a parametric study for the energy requirements [108].They reported both individual and multiple effects of parameters such as MVD, LWC, surface temperature, ambient temperature, airspeed, and angle of attack on the energy requirement.Zhou et al. investigated the effect of the meteorological conditions on the surface temperature and the runback ice [113].They concluded that the LWC and the Mach number had more influence than the air temperature.Papadakis et al. investigated the effect of piccolo parameters, including its geometry, hot-air temperature, and hot-air mass flow rate, on the hot-air system performance.They confirmed that the surface temperature was sensitive to the geometric relation between the piccolo tube, the leading edge of the wing, and the hole pattern on the piccolo tube [114].They also reported the effect of the hot-air mass flow and temperature on the surface temperature [115].Zhang et al. conducted the sensitivity analysis on the hot-air system performance [116].They concluded that the hole diameter on the piccolo tube was one of the critical structural parameters.In addition to optimizing thermal icing mitigation for the wing, a similar trend can be seen to optimize the mitigation of icing on the engine inlet guide vane [117][118][119][120].
(2) Mechanical Method The mechanical methods refer to the mitigation techniques that apply a mechanical force onto the accreted ice to break and remove ice from the surface.Since it allows a certain amount of ice accretion, this icing mitigation technique is considered deicing.One of the representative mechanical mitigation systems is a pneumatic boot (Figure 9).A sheet of rubber is attached to the surface of the leading edge of the wing, and it breaks the accreted ice when inflated by pressurized air.Compared to the thermal method, this system has advantages in small weight and low cost [121].Therefore, aircraft such as helicopters and propeller airplanes adopt pneumatic boots due to the limited weight tolerance and energy source.Since the pneumatic boot deicer can remove the accreted ice only after the thickness becomes larger than a threshold value, an improvement in its efficiency is required so that the ice can be removed at as thin a threshold as possible [122].As more efficient deicing systems were demanded, several electromechanical methods were proposed [122].The electroexpulsive separation system (EESS) consists of pairs of two layers of boots (Figure 10).When current flows in opposite directions between the layers, the magnetic fields produce force such that the layers separate from each other, resulting in the deformation of the outer layer.Accreted ice breaks up due to this deformation.The electromagnetic impulse deicer (EIDI) uses coils embedded under the metal skin and a high-voltage capacitor (Figure 11) [123].When the capacitor is discharged, the coils create a rapidly forming and collapsing electromagnetic field.The magnetic field induces the eddy currents in the metal skin, causing a repulsive force that detaches the ice [124].(3) Chemical Method A mitigation technique using a chemical method has been studied.It utilizes a substance's chemical properties to either lower freezing temperature or lower ice adhesion to prevent icing.The former approach can be seen on the ground before the flight [126].Depending on the mixture of chemical fluids, this approach works as both deicing and anti-icing.The freezing depressant for deicing consists of ethylene or propylene glycol and hot water and is used to remove the accreted ice [127].On the other hand, the fluid for anti-icing is similar to that for deicing, and it is known that polymers increase the viscosity [127].
While the chemical method has been used for ground de-/anti-icing, the environmental influence of those freezing depressants has been an area of concern [128,129].When a glycol product is released to the soil or river near airports, the biochemical oxygen demand (BOD) increases.Here, BOD is the amount of oxygen consumed by bacteria or microorganisms to decompose organic materials and is used as one indicator for the quality of water.Several studies on the environmental effect of freezing depressants have been conducted, and it was revealed that additives in the freezing depressant rather than glycol contributed to the increase in BOD [130,131], especially benzotriazole, which inhibits metal corrosion [132].Since then, investigations of the influence of freezing depressants on the environment have focused on the presence of benzotriazole [133][134][135][136].

Passive Mitigation
(1) Physicochemical Method The next set of techniques used to mitigate ice adhesion uses a physicochemical approach.It reduces ice adhesion on an aircraft surface by changing the surface properties so that ice can be easily removed from the surface [137].This surface modification has been mainly done through surface coatings, and those coatings are often called icephobic coatings.Although the development of icephobic coating has been conducted over decades [138], they have not yet been adopted by commercial aircraft.The focus of the research is the low ice adhesion as well as the durability of the icephobic surface [139].
Hydrophobic coatings have long been investigated for their icephobicity [140].On hydrophobic coatings, ice formation is delayed since water droplets are repelled before crystallization [141].In addition, the smaller contact area between the water droplet and surface results in a low ice adhesion [140].Developments of hydrophobic coatings as icephobic coatings are still ongoing [142][143][144][145][146][147][148][149].Recently, the low ice adhesion of 9 kPa was reported on a micro-nanostructured surface in [150].Since the 2010s, the slippery lubricantinfused porous surface (SLIPS) has gained attention as an icephobic coating [151,152].This type of coating works as a lubricant between the ice and surface, resulting in the low ice adhesion between 15 and 25 kPa reported in [153,154] compared to the ice adhesion of 821.9 kPa for bare aluminum [155].However, the icephobicity of SLIPS had a short lifetime since the lubricant was depleted through evaporation or consumption [156].In order to increase the durability, an icephobic surface with bioinspired solid organogel materials has been developed [157].While the sacrificial alkane surface layer is depleted with the removed ice, the layer can be regenerated from the crosslinked poly-dimethylsiloxane (PDMS) swelled by the molten alkane.More details on icephobic coatings can be found in the summaries of current studies presented by Shen et al. [158] and Zhuo et al. [159].

Hybrid Mitigation
Icephobic coatings are an ideal mitigation technique in energy consumption since they do not require any energy source.While various icephobic coatings have been developed, none of the coatings has been proven applicable for aerospace applications [160].Instead, the hybrid mitigation technique in which the coatings are used to support an active mitigation technique has been considered.One such example is the mitigation technique combining electrical heating and an icephobic coating proposed by Morita et al. [161].To lower the power consumption of the electrical heater, the heating temperature should be lower and/or the heating area should be minimized.However, in that case, ice accretion occurs when impinged water is moved to the unheated region by the local airflow.The icephobic coating based on a hydrophobic surface is applied right after the heating area, which prompts the detachment of the water from the surface of the airfoil before freezing.They reported that the combination of the thermal method and the icephobic coating only used 30-70% of the power consumption of the thermal method itself.Another hybrid mitigation was proposed in [162] which combined the icephobic coating and electromechanical mitigation techniques such as EESS and EIDI.Although the concept has been proposed, further investigation is required to evaluate its effectiveness due to the low maturity stage of the research [160].

Conclusions
The icing on an aircraft has been a major threat to flight safety.Researchers have focused on understanding the underlying physics of icing and developing better icing mitigation techniques in energy consumption.This review article aimed to summarize such advancements in icing physics and icing mitigation.Research on the icing physics was reviewed in three parts: fixed wings, rotors, and engines.Researchers focused on microscale phenomena such as droplet splashing and surface roughness of the accreted ice for icing on a fixed wing.For icing on rotors, implementing three-dimensional effects such as the ice shedding into icing simulations has been the main focus of the research.Engine icing, also known as the ice crystal icing, was a relatively new area of the research compared to icing on fixed wings and rotors.Researchers worked to understand the fundamental mechanism of ice crystal icing and develop better simulation tools.
Five types of icing mitigation methods were reviewed as the current technology in aircraft icing.The thermal and mechanical methods as active mitigation techniques have already been used in flight operations.For the thermal method, there has been a significant push to achieve lower energy consumption.On the other hand, several new types of mitigation techniques using electrical power were proposed for the mechanical method.The chemical method as passive mitigation has been used to mitigate the icing on the aircraft on the ground.Since the freezing depressant may influence the environment, considerable research activities have been carried out to evaluate the environmental impacts of this ice mitigation technique.An additional method, the physiochemical method, was developed as a passive mitigation technique.Hydrophobic coatings have been developed as icephobic coatings, and recent interest has grown surrounding the slippery lubricant-infused porous surface (SLIPS).The hybrid mitigation technique was also introduced.It lowers the energy consumption of the thermal method by combining it with the icephobic coating.
As can be seen from this review of the current achievements in aircraft icing physics and techniques for the mitigation of aircraft icing, the focus of most research has been on (1) acquiring the knowledge on details of underlying icing physics and (2) improving the mitigation techniques in terms of energy consumption.However, the ice crystal icing knowledge seems immature due to its short history of research activity.Consequently, the primary mitigation of the ice crystal icing is to avoid the area where the ice crystal icing is likely to occur.More investigations are demanded, especially for the ice crystal icing, to achieve better airworthiness.For an ice mitigation technique, hybrid mitigation is introduced, combining an active and a passive technique.It demonstrated a reduction in energy consumption compared to that of the existing active mitigation.Based on the current review, the research on this mitigation technique may grow.New hybrid mitigations will be seen, and their research challenges need to be highlighted.

Figure 1 .
Figure 1.Schematic of parameters for the mass conservation in Messinger model.

Figure 2 .
Figure 2. Main splashing parameters related to the numerical models.
29 h of testing time per day on average.Due to the high cost and safety concerns for in-flight icing tests, rotor icing experimental facilities on the ground were established in several regions.The National Research Council Canada constructed a spray rig outside, producing water droplets between 20 and 60 µm in diameter.The U.S. Air Force has McKinley Climatic Chamber in Florida [66].The Pennsylvania State University has a smaller test stand in the Adverse Environment Rotor Test Stand (AERTS) where rotors up to 4.5 ft in diameter fit in the test section (Figure 3) [67].

Figure 5 .
Figure 5. Structure of the engine and potential ice accretion area (reproduced from [7]).

Figure 7 .
Figure 7. Schematic of parameters contributing to the mass conservation (obtained from [96]).Note that the impinging mass rate, .m imp , can be used to describe the impinging water