# Modeling in-flight ice accretion under uncertain conditions

Giulio Gori, Pietro Marco Congedo, Olivier Le Maitre, Tommaso Bellosta,
Alberto Guardone

► To cite this version:

Giulio Gori, Pietro Marco Congedo, Olivier Le Maitre, Tommaso Bellosta, Alberto Guardone. Modeling in-flight ice accretion under uncertain conditions. Journal of Aircraft, 2022, 59 (3), 10.2514/1.C036545. hal-03707307v2

HAL Id: hal-03707307

https://inria.hal.science/hal-03707307v2

Submitted on 28 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling in-flight ice accretion under uncertain conditions

G. Gori, P.M. Congedo

INRIA, CMAP, École Polytechnique, IPP, Route de Saclay, 91128 Palaiseau, France.

O. Le Maître

INRIA, CNRS, CMAP, École Polytechnique, IPP, Route de Saclay, 91128 Palaiseau, France.

T. Bellosta, A. Guardone

Department of Aerospace Science & Technology, Politecnico di Milano, Via La Masa 35, 20156 Milano, Italy

In-flight ice accretion under parametric uncertainty is investigated. Three test cases are presented which reproduce experiments carried out at the NASA's Glenn Icing Research Tunnel (IRT) facility. A preliminary accuracy assessment, achieved comparing numerical predictions against experimental observations, confirm the robustness and the predictiveness of the computerized icing model. Besides, sensitivity analyses highlight the variance of the targeted outputs with respect to the different uncertain inputs. In rime icing conditions, a predominant role is played by the uncertainty affecting the airfoil angle of attack, the cloud liquid water content and the droplets' mean volume diameter. In glaze icing condition, the sensitivity analysis shows instead that the output variability is due mainly to the ambient temperature uncertainty. Results expose a major criticality of standard uncertainty quantification techniques. The issue is inherent in the approximation of the full icing model behavior in domain regions scarcely affected by ice build up. To mitigate the issue, a non-linear regression method is proposed and applied.

## Nomenclature

$A$ = elemental cell surface area [m²]

$AoA$ = angle of attack [deg]

$B$ = ice thickness [m]

$B_g$ = rime ice limit thickness [m]

$C_D$ = drag coefficient [-]

$C_p$ = specific heat at constant pressure [J/kgK]

$d$ = diameter [m]

$E^t$ = specific total energy [J·m³/Kg]

$g$ = gravitational acceleration [m/s²]

$IRT$ = icing research tunnel
$k$ = thermal conductivity [W/s²]

$LWC$ = liquid water content [g/m³]

$L_F$ = ice latent heat of fusion [J/Kg]

$L_S$ = ice latent heat of sublimation [J/Kg]

$MVD$ = mean volume diameter [µ·m]

$m$ = mass [Kg]

$\dot{m}_d$ = mass flux due to impinging droplets [Kg/s]

$\dot{m}_{in}$ = mass flux due to liquid water incoming from neighbouring cells [Kg/s]

$\dot{m}_{out}$ = mass flux due to liquid water leaving to neighbouring cells [Kg/s]

$PCE$ = Polynomial Chaos Expansion

$QoI$ = Quantity of Interest

$\dot{Q}_{bottom}$ = Heat flux exchanged between the phase changing interface and ice [W/m²]

$\dot{Q}_{top}$ = Heat flux due to sublimation [W/m²]

$\dot{Q}_{top}$ = Heat flux exchanged between the phase changing interface and water (or air) [W/m²]

$Re$ = Reynolds number [-]

$s$ = curvilinear abscissa [m]

$T_\infty$ = freestream temperature [m/s]

$t$ = time [s]

$U$ = vector of conserved variables

$V_\infty$ = freestream velocity [m/s]

$v$ = local velocity vector [m/s]

$We$ = Weber number [-]

Greek letters

$\beta$ = collection efficiency [-]

$\mu$ = viscosity [Pa·s]

$\rho$ = density [Kg/m³]

$\tau$ = viscous stress tensor

Subscripts

air = air

d = droplet

ig = glaze ice

ir = rime ice
I. Introduction

In-flight ice accretion is a recurring detrimental phenomenon leading to a modification of the airflow around wings, flight sensors and engines, and to an increase of the overall aircraft weight. Ultimately, it results in performance degradation and poses a serious threat to safety. Starting from the pioneering work of J. Stefan [1], originally devoted to the study of ice formation in the arctic ocean, the last decades have witnessed the development of a large number of mathematical models devised for investigating in-flight icing, also due to the high cost of ice wind-tunnel experiments and in-flight testing. In 1953, Messinger proposed the first accretion model specifically tailored to aerospace applications [2]. Since then, growing concerns from regulators and manufacturers to ice-related safety issues fostered the related research activity. Practically, the icing modeling activity always resulted into the implementation of a wide variety of numerical codes, see e.g. [3–15] Computational models help investigating ice build-up, providing a mean to predict ice formation at critical locations over the aircraft fuselage and wings. Ultimately, codes are exploited to complement the experimental activity required to design and certify new aircraft [16].

Unfortunately, our understanding of reality is unavoidably affected by uncertainty which questions the credibility of computational models and, ultimately, the reliability of numerical predictions. This is particularly true for in-flight ice accretion, which abides by complex physical laws challenging standard Uncertainty Quantification (UQ) techniques. Being capable of generating predictions under uncertainty is of the utmost relevance for enhancing the robustness of in-flight ice protection systems and the aircraft certification process. At the time of the writing, published research concerning UQ in icing applications is seldom. In [17], a standard Polynomial Chaos (PC) approach and a standard Monte Carlo method (MC) were employed to investigate the effects of different ice structures w.r.t. the aerodynamic performance of 2D airfoils. By means of a parametric description of some canonical ice shapes, the UQ analysis targets the aerodynamic performances rather than the accretion mechanism itself or the shape of the ice crystal forming over the wing. Recently, authors [18] started evaluating the effects of some input uncertainty e.g., ice surface roughness and density, w.r.t. specific features characterizing the ice shape resulting from the accretion process i.e., ice accretion limits, horns position, horns angle and others. The authors take advantage of a standard Radial Basis Function (RBF) approach to build a surrogate of the high-fidelity computational model, to ultimately evaluate the output sensitivity via the Sobol’s sequence sampling method.

The present paper aims at quantifying the sensitivity of the predicted ice shape to uncertain test conditions with a twofold objective. On one hand, the aim is to shed light on the physics by exposing the relative importance of the diverse uncertain inputs to the resulting ice shape. On the other hand, this paper is also meant to be an embryonic effort towards assessing the accuracy of numerical predictions against real observations. Three different test cases, concerning
<figure><img src="image_2.png" /><figcaption>Fig. 1 The flowchart of the computational model.</figcaption></figure>

experiments carried out in the NASA Icing Research Tunnel (IRT) facility, are reproduced numerically, using the PoliMIce ice accretion simulation framework [6]. The first test case targets the distribution of water collected over the wing of an business aircraft flying across a cloud [19]. The remaining two test cases target instead the ice layer thickness and aim at reproducing the experiments described in [20]. In a first stage, the input uncertainties are forward-propagated across the numerical framework, to selected Quantities of Interest (QoI), using a standard (linear) Polynomial Chaos expansion approach [21]. Results from standard PC expose a major criticality inherent the approximation of the full icing model at locations scarcely subject to ice accretion. To mitigate the issue, we propose and apply a non-linear PC approach to approximate continuous functionals achieving a constant lower bound on some part of their support, comparing the resulting approximation against the one obtained by using standard techniques.

The paper is structured as follows. In Section II, the computational model employed to simulate ice accretion is presented. Section III briefly presents the classical PC approach to uncertainty quantification and the non-linear regression method we propose. In Section IV, the input uncertainties are propagated to the model outputs and numerical predictions are compared against experimental data. Section V reports the final comments and hints to possible future research.

## II. The computational model for in-flight ice accretion

In this Section we present the computational model we employ to predict ice accretion in aeronautical applications. It consists in a multi-disciplinary framework loosely coupling an in-house Lagrangian particle tracking code with the open-source Computational Fluid Dynamics (CFD) solver SU2 [22] and the ice accretion solver PoliMIce [6]. The framework flowchart is reported in Fig.1.

The computational model approximates clouds as uniform, homogeneous and isotropic. The problem is treated as a multi-phase dispersed mist flow around the airplane. With reference to Fig. 1, first the aerodynamic field developing around the airfoil, or more generally the aircraft, is reconstructed using the SU2 RANS solver described in Sec. II.A.
Then, the field solution is passed to our in-house particle-tracking software, PoliDrop, see Sec. II.B. The solver computes the trajectory of each single droplet, from its initial position within the cloud at the far-field ahead of the aircraft to the impinging point over the wings and fuselage, in order to estimate the amount of water collected by the aircraft surface. In other words, the particle tracker is exploited to compute the so-called collection efficiency ($\beta$). The ice thickness is then computed from the distribution of $\beta$ and additional aerodynamic quantities including e.g. the aircraft velocity, the temperature and the shear stresses, using a liquid-solid phase transition model which accounts also for additional uncertain parameters (the ice thermal conductivity, the ice density and so on), see Sec.II.C. This transition model is based on an equilibrium energy balance and it is implemented within the PoliMIce ice accretion engine.

## A. The aerodynamic field model

The aerodynamic field developing around the geometry under investigation is reconstructed by means of the well-known Reynolds-Averaged Navier-Stokes (RANS) model for compressible fluid flows. To the purposes of this paper, the computational domain consists of a two-dimensional airfoil immersed in a free stream flow, see Fig. 2 reporting the Mach number field developing around the targeted geometry. We consider a domain $\Omega \subset \mathbb{R}^2$ with a

<figure><img src="image_3.png" /><figcaption>Fig. 2 The computational domain employed to reconstruct the aerodynamic field surrounding the airfoil. The picture reports the Mach number field.</figcaption></figure>

disconnected boundary divided into $\Gamma_{\infty}$, the far-field component, and $\Gamma_S$, an adiabatic wall boundary condition applied on the airfoil. The differential form of the RANS equations reads

$$ egin{cases} 
\mathcal{R}(\mathbf{U}) = \frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}^c - \nabla \cdot (\mu^{\nu k} \mathbf{F}^{\nu k}) - (Q) = 0 & \text{in } \Omega, \quad t > 0, \\ 
\mathbf{v} = \mathbf{0} & \text{on } \Gamma_S, \\ 
\partial_n T = 0, & \text{on } \Gamma_S, \\ 
(\mathbf{W})_+ = \mathbf{W}_{\infty} & \text{on } \Gamma_{\infty}, 
\end{cases} 
\quad (1) $$
being $\mathbf{U} = \{\rho, \rho\mathbf{v}, \rho E^t\}^T$ the vector of the conserved variables whereas $\mathbf{F}^c, \mathbf{F}^{vk}$ and $\mathbf{Q}$ are, respectively, the convective and the viscous fluxes and the source term, reading

$$ \mathbf{F}^c = \left\{ \begin{matrix} \rho \mathbf{v} \\ \rho \mathbf{v} \otimes \mathbf{v} + \frac{\mathbf{I} \mathbf{P}}{\mathbf{\Xi}} \\ \rho E^t \mathbf{v} + \mathbf{v} \end{matrix} \right\}, \quad \mathbf{F}^{v1} = \left\{ \begin{matrix} \mathbf{.} \\ \mathbf{\tau} \\ \mathbf{\tau} \cdot \mathbf{v} \end{matrix} \right\}, \quad \mathbf{F}^{v2} = \left\{ \begin{matrix} \mathbf{.} \\ c_p \nabla T \\ c_p \nabla T \end{matrix} \right\}, \quad \mathbf{Q} = \left\{ \begin{matrix} \mathbf{.} \\ \mathbf{\tau} \\ \mathbf{\tau} \cdot \mathbf{v} \end{matrix} \right\}. \qquad (2) $$

with $\mu$ dynamic viscosity coefficient and $\mathbf{\tau}$ viscous stress tensor. Note that in the above expressions we dropped the subscript 'air' to lighten the notation, but naturally all the quantities refer to the carrier fluid which is standard air.

The aerodynamic field depends on a set of parameters, imposed by means of the term $\mathbf{W}_\infty$ in Sys. (1), describing the environment (freestream air pressure $P_\infty$ and freestream air temperature $T_\infty$) and specifying the aircraft flying conditions (freestream Mach number $M_\infty$ and Angle of Attack AoA). These parameters vary significantly as the flight proceeds through the planned route, crossing diverse geographical regions and undergoing mass variations (which requires adjustments of the aircraft attitude) due to fuel consumption. Nevertheless, such variations generally occur over a large time interval, surely larger than the ice accretion time scale, therefore justifying the assumption of an averagely steady airflow. Unsteady turbulent effects are accounted for according to the Boussinesq hypothesis. The turbulent viscosity is obtained via the Menter's shear-stress transport model [23].

Within the SU2 solver, the compressible Reynolds averaged Navier-Stokes equations are discretized using the finite volume method with a standard edge-based structure on a dual grid with control volumes constructed using a median-dual vertex-based scheme. Convective fluxes are discretized using a limited second order MUSCL scheme with an Approximate Riemann Solver (ARS) of Roe type. In our setting, the Venkatakrishnan flux limiter is employed. Viscous fluxes are discretized using a standard average of gradients approach. Source terms are approximated at each node using a piecewise-constant reconstruction within each control volume. Gradients are obtained via a weighted least-squares approach. A time-marching approach is used to drive the RANS system to a steady solution using an implicit Euler scheme. A multi-grid approach is employed to accelerate the convergence of the CFD solver iterations to the steady solution. Convergence is monitored by checking the root mean square of the density equation residual. Computations are stopped once the residual is at least eighth orders of magnitude smaller than the first iteration.

## B. The cloud model

A fundamental parameter characterizing clouds is the so-called Liquid Water Content (LWC), expressing the spatial concentration of droplets suspended in the atmosphere. Knowing the concentration of dispersed particles is key for determining the coupling between phases in particle laden flows. The LWC measures the amount of water (mass) specified to the unit volume of air, and it varies significantly depending on the type of cloud encountered, with typical
values ranging from 0.002 to 3 g/m³. Given the low water volume fractions of typical clouds, a one-way coupling can be generally assumed for the particle laden flow involved in-flight ice-accretion problems. Therefore, the carrier flow equations are solved independently from the presence of a dispersed phase. Once the carrier field is available, see Sec. II.A, the in-house Lagrangian solver PoliDrop is exploited to reconstruct droplet trajectories and, ultimately, to determine the impinging points over the aircraft surface. The dispersed phase is handled by means of the Discrete Parcel Method (DPM) [24]: each particle, or parcel, represents a set of neighboring droplets and it is characterized by an average velocity $\mathbf{u}_p$ and an average diameter $d_p$. In the Lagrangian setting, the trajectory of each parcel is reconstructed by numerically integrating the following differential equation

$$m_p \frac{d\mathbf{v}_p}{dt} = \frac{\pi}{8} \mu_{\text{air}} d_p \text{Re}_p (\mathbf{v}_{\text{air}} - \mathbf{v}_p) C_D + \frac{\pi}{6} d_p^3 \mathbf{g} (\rho_p - \rho_{\text{air}}), \quad (3)$$

where $\mu_{\text{air}}$, $\rho_{\text{air}}$, $\mathbf{v}_{\text{air}}$ are the air local viscosity, density and velocity, $C_D$ is the parcel drag coefficient and $\mathbf{g}$ is the gravitational acceleration. Re$_p$ is the Reynolds number computed using the relative droplet-flow velocity

$$\text{Re}_p = \frac{\rho_{\text{air}} (\mathbf{v}_{\text{air}} - \mathbf{v}_p) d_p}{\mu_{\text{air}}}. \quad (4)$$

Note that, since water has a much larger density than air ($\rho_p \gg \rho_{\text{air}}$) [25], the only forces considered in the equation of motion (3) are gravity and aerodynamic drag. The drag coefficient model employed in our setting is taken from Morrison [26]. This model takes into account deformations deviating the particle shape from a perfect sphere. Such feature is fundamental to simulate the motion of droplets of large dimensions, For these latter, the deformation can be approximated to different degrees in between the shape of a perfect sphere and an oblate disk. The parameter quantifying the extent of this eccentricity is defined as $f = (1 + 0.07\sqrt{\text{We}})^{-6}$, see [27], being We the Weber number

$$\text{We} = \frac{\rho_p d_p (\mathbf{v}_p - \mathbf{v}_{\text{air}})^2}{\sigma_p}, \quad (5)$$

relating the drag force acting on the droplet to its surface tension $\sigma_p$. The drag coefficient of the deformed droplet can be obtained by means of a weighted average of that of a rigid sphere and that of a disk:

$$C_D = \begin{cases} (1 - f)C_{D\text{Sphere}} + fC_{D\text{Disk}} & \text{if } \text{We} \le 12 \\ C_{D\text{Disk}} & \text{if } \text{We} > 12 \end{cases} \quad (6)$$

where the disk drag coefficient is reported in [28].

In our model, clouds are represented as a single front initially placed at an arbitrary distance ahead of the aircraft. This distance is set so that droplets are traced starting from an unperturbed region of the domain and so that the
computational burden related to the trajectory integration in time, proportional to the distance from the surface, is
limited. In three dimensional problems, this front consists of a two-dimensional layer of droplets uniformly and
randomly distributed over a plane oriented perpendicularly to the free stream direction. In a two-dimensional setting,
which is the case considered hereinafter, the layer degenerates into a straight line normal to the mean flow direction. In
practical applications, clouds are made up by a massive amount of droplets of diverse size. We assume here that the
front only contains droplets of constant size, set equal to the cloud Median Volume Diameter (MVD) i.e., the diameter
corresponding to the median value of the distribution characterizing the size of droplets contained in the cloud.

As mentioned, by integrating Eq. (3) numerically it is possible to reconstructs the trajectory of each droplet from
its initial position, at the cloud front, to its potential impinging point over the aircraft surface. In aeronautical icing
applications, the impingement of a droplet against a solid surface can result in four different outcomes [29] namely,
stick contact (complete adhesion and deposition), rebound (complete bouncing), spread (on the surface) and splash
(partial deposition with formation of secondary droplets). The modeling of the droplet-wall interaction upon impact
is of the paramount importance as it allows to accurately estimate the amount of water collected over the surface, to
ultimately evaluate the thickness of the resulting ice layer. The occurrence of a particular mechanisms depends on many
parameters such as the droplets characteristics (size, density, viscosity, surface tension), the impact velocity and the
impact angle, and the properties of the surface. For a thorough description of the impingement models used, the reader
is referred to Ref. [30].

In aeronautical icing problems, the collection efficiency parameter $\beta$ is typically employed to express the fraction of
water mass, contained in the free stream, collected at a given location over the aircraft surface. It can be interpreted as
the ratio of the cross section area of a droplet stream tube at the farfield to the area of impact on the body surface, as
reported in Fig. 3. Alternatively, it can be seen as the ratio of the surface collection density i.e., the number of droplets
collected, at a given location over the aircraft per unit surface (or per unit length), w.r.t. the nominal cloud front density.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Collection efficiency parameter $\beta$ computed for a two and a three-dimensional problem, from Ref. [16]</figcaption></figure>
### C. The phase-change model

Depending on the flight conditions, ice can accrete on the aircraft surface according to different mechanisms leading to the formation of characteristic structures and layer shapes made up of ice with different physical properties. Typically, the ice accretion regime is classified as either *rime* or *glaze*. Rime ice occurs at very low temperatures and it results from an instantaneous liquid-to-solid phase transition. Namely, at very low temperatures clouds may be composed by super-cooled droplets i.e., small liquid water spheres characterized by a temperature lower than the freezing point. This peculiar counter intuitive state corresponds to an unstable thermodynamic equilibrium point. When the aircraft flies across the cloud, droplets hit the surface and their unstable equilibrium state is perturbed. Therefore, an instantaneous water-to-ice phase transition is triggered. The particularly short time scale upon which the phase-change occurs causes small air bubbles to remain trapped within the ice layer. As a result, the ice is characterized by an opaque aspect and a low density. Differently, glaze ice forms at higher temperatures following a slower freezing process. At such conditions, droplets impact the surface and remain there forming a thin liquid film before freezing. These particular circumstances allows for air bubbles to separate from water, giving ice an homogeneous denser structure which results into a typical transparent and lucid look. Ice accretion of mixed type can also possibly occur, with both rime and glaze ice accreting one after the other.

From a practical point of view, computing the thickness of the ice layer implies solving a phase-change problem over the domain solid boundaries i.e., the aircraft surface. Typically, boundaries are first discretized into a set of small elemental cells of simple shape. Within each sub-domain, a one-dimensional phase-transition problem is solved (in the direction of the local surface normal). Figure 4 illustrates an exemplary elemental cell. Sub-domains are loosely coupled through boundary conditions enforcing the conservation of mass at cells interfaces. Namely, through specific terms representing the amount of unfrozen water entering/leaving the cell ($\dot{m}_{in}/\dot{m}_{out}$) from/to the upstream/downstream elemental domain, under the action of the external flow force. Ice grows in the direction normal to the wall, due to a certain amount of freezing water collected from the cloud ($\dot{m}_d$) and from the adjacent cell ($\dot{m}_{in}$). The terms $\dot{Q}_{top}$ and $\dot{Q}_{down}$ represents instead thermal power contributions. Namely, $\dot{Q}_{top}$ is a generic term including for instance the amount of power dissipated to the external airflow through convection, the amount of power absorbed due to the droplets kinetic energy release or due to the air friction and many other contributions. Instead, $\dot{Q}_{down}$ practically accounts for the thermal power dissipated via conduction thorough the ice layer, to the aircraft structure.

In this work, to simulate ice accretion we employ the PoliMIce accretion model implementing a local solution to the unsteady Stefan problem [7]. It consists in a modified version of the Myers' extension to the Messinger's model, explicitly accounting for the mass fluxes related to sublimation and to a more detailed description of the liquid film flow above the ice surface. Moreover, it includes an unsteady description of the heat diffusion problem within the ice layer based on the exact local solution of the unsteady Stefan problem for the temperature profiles within the ice layer in glaze conditions, see [7]. Summarizing, the icing model employed in this work recognizes the rime and glaze regimes,
<figure><img src="image_5.png" /><figcaption>Fig. 4 Elemental domain for the computation of the ice layer thickness over the surface of interest.</figcaption></figure>

providing the following accretion laws:

$$ \text{Rime:} \quad \frac{\partial B}{\partial t} = \left[ \frac{\dot{m}_d(\beta) + \dot{m}_{in} - \dot{m}_{out}}{A \rho_{ir}} \right], \quad \text{if } B < B_g, \quad (7) $$

$$ \text{Glaze:} \quad \frac{\partial B}{\partial t} = \frac{1}{\rho_{ig} L_F} \left[ k_{ig} \frac{(273.15 - T_{\text{Wall}}) e^{-\lambda^2} \sqrt{\rho_{ig} C_P}}{\text{erf}(\lambda) \sqrt{\pi k_{ig} t}} + \frac{\dot{Q}_{\text{top}}(\beta)}{A} \right], \quad \text{if } B \ge B_g. \quad (8) $$

where $\lambda(t) = \frac{B(t)\sqrt{\rho_{ig}C_P}}{\sqrt{k_{ig}}}$. Note that, according to standard practice, we assume $T_{Wall}$ to be equal to the ambient temperature. The rime-to-glaze accretion transition is triggered by the reaching of a particular ice layer thickness acting as a threshold separating the two accretion regimes. This limiting thickness is computed as

$$ B_g = \frac{Ak_{ir}(273.15 - T_{\text{wall}})}{2 \left[ L_F (\beta \text{LWC } V_{\infty} A - \dot{Q}_s L_s^{-1}) - \dot{Q}_{\text{top}}(\beta) \right]}. \quad (9) $$

Whenever the thickness of the ice layer is lower than $B_g$, accretion occurs in the rime regime. Otherwise, glaze ice appears. The reader is referred to the nomenclature section for the definition of each entry included in the aforementioned equations and to Ref. [7] for a comprehensive analysis of the accretion model.

All the equations described above depend implicitly on the parameter $\beta$, which varies in time as the ice shape changes giving birth to ridges and valleys across the frozen water layer, ultimately affecting the aerodynamic field in the close proximity of the aircraft surface. Note that, to the purposes of this paper, numerical simulations are obtained using the so-called one-shot setting. That is, the aerodynamic field and the collection efficiency are computed once, at the initial time instant and considering the dry profile. After, the ice accretion engine computes the ice layer thickness forming after a given interval considering a single time step. Therefore, the collection efficiency is assumed to be a steady parameter. Nevertheless, we point out that, even if multi-step simulations are not employed in this work, they are possible within our computational framework. In addition to the $\beta$ parameter, ice accretion also depends on the value of LWC, on the ice density $\rho_{ir}$ and thermal conductivity $k_{ir}$ and the ice latent heat of fusion $L_F$. These latter are physical
properties of ice and, as such, vary depending on the structure of the crystals forming as a consequence of the freezing process. Due to the possible presence of small air bubbles trapped within the ice layer, significant differences are found in between the structure of crystals of the rime and glaze type. Empirical correlations exist to compute the values of both $\rho_{ir}$ and $k_{ir}$, depending on the accretion regime. The Jones' model [31] can be used to estimate the ice density starting from the problem conditions and accounting for the geometrical characteristic of the airfoil (which affect the droplet-wall interaction upon impact). The density is computed as:

$$ \rho_{ir} = 1000 \exp(-0.15(1 + 6043S^{-2.65})), \quad \text{with} \quad S = -\frac{\text{MVD}^{0.82}\text{V}_{\infty}^{0.59}\text{LWC}^{0.21}}{\text{D}^{0.48}(T_{\infty})^{0.23}}, \qquad (10) $$

where D is the mean cylinder diameter resulting from the rotating multi-cylinder method used by Jones to obtain the data, see [31]. In aeronautical applications, the definition of a mean diameter for an airfoil is ambiguous. Therefore, the surface curvature at the leading edge is used instead. The dependency of the ice thermal conductivity on density can again be modeled using experimental correlations [32]

$$ k_{ir} = [2.107 + 0.003618(\rho_{air} - \rho_{ir})] \hat{k}_{i} / \hat{k}_{i}^{\text{ref}} \qquad (11) $$

where a multiplicative factor $\hat{k}_i/\hat{k}_i^{\text{ref}}$ is applied to recover the dependence with temperature. Indeed, the correlation expression (11) was established for an ambient temperature of -3 Celsius degrees. In the multiplicative factor $\hat{k}_i/\hat{k}_i^{\text{ref}}$, $\hat{k}_i$ is the thermal conductivity of pure ice at the target temperature whereas $\hat{k}_i^{\text{ref}}$ is that of pure ice at the reference temperature for which the formula was devised (-3° Celsius). According to [32], the above formula yields accurate results for a porous ice density between 550 and 917 kg/m³, values typical of in-flight ice applications.

Clearly, in view of Eq. (7), Eq. (8) and Eq. (9), ice accretion is a highly non-linear and complex mechanism. The rime-to-glaze transition is triggered according to a threshold function introducing discontinuities which contribute in increasing the complexity of the model response to different inputs and conditions. This complex mechanism often results in a lack of smoothness of the targeted QoI with respect to the input parameters.

### III. Uncertainty Quantification

To forward-propagate the uncertainty affecting test conditions (x) to the targeted QoIs (y = f(x), a generic output of the computational model), we employ a Polynomial Chaos (PC) expansion approach [21, 33]. The PC expansion consists in a surrogate of the mapping x $\mapsto$ f(x) resulting from the evaluation of the full computational model. By assuming that the uncertain inputs are parametrized using a finite set of d independent random variables $\xi = \xi(x) = (\xi(x)_1, \dots, \xi(x)_d) \in \Omega^d$ with known probability densities p(\xi), the expansion corresponds to the orthogonal
projection, into the linear span, of a finite set of orthonormal polynomials $\Psi_k(\xi)$

$$f(\xi) = \sum_{k=0}^{\infty} c_k \Psi_k(\xi), \qquad (12)$$

where $c_k$ are the PC coefficients and where the index of summation $k$ indicates the polynomial order. $k$ ranges from 0 to an upper arbitrary value $P$ set to truncate the expansion. The polynomial basis $\Psi(\xi)$ is chosen such that the following orthogonality relation, corresponding to the inner product in the $L_2$ Hilbert space, is satisfied

$$\langle \Psi_i(\xi), \Psi_j(\xi) \rangle = \int \Psi_i(\xi) \Psi_j(\xi) p(\xi) d\xi = \langle \Psi_i(\xi)^2 \rangle \delta_{ij}. \qquad (13)$$

Thanks to properties of PC, the uncertain model output can be characterized by manipulating the coefficients of the truncated expansion. The output mean $\bar{y}$ and variance $\sigma^2$ can indeed be approximated as

$$\bar{y} = \mathbb{E}[y] \approx c_0, \qquad \sigma^2 = \mathbb{E}[(y - \bar{y})^2] \approx \sum_{k=1}^{P} c_k^2. \qquad (14)$$

Similarly, Sobol indexes for sensitivity analysis can be straightforwardly obtained since the functional dependencies, w.r.t. the input, are explicit, see [34]. The Sobol functional decomposition [35] relies on the assumption of that

$$f(\xi) = \sum_{u \le \{1, 2, \dots, d\}} f_u(\xi_u), \qquad (15)$$

with $s = \text{card}(u) = |u|$ and $f_0 = f_0$. Each $f_u$, except $f_0$, holds to the relation

$$\int f_u(\xi_u) p(\xi_u) d\xi_u = 0, \quad \forall u \ni i. \qquad (16)$$

By defining the output variance $\sigma^2$ and the conditional variances $\sigma_u^2$ i.e., the variances resulting from the Sobol decomposition,

$$\sigma^2 = \int_{\Omega^d} f^2(\xi) p(\xi) d\xi - f_0^2(\xi), \qquad \sigma_u^2 = \int_{\Omega^{|u|}} f_u^2(\xi_u) p(\xi_u) d\xi_u, \qquad (17)$$

it is possible to define the so-called Sobol sensitivity indexes

$$S_u = \frac{\sigma^2}{\sigma_u^2}, \quad \text{with} \quad \sum_{u=\{1, 2, \dots, d\}} S_u = 1 \quad \text{and} \quad u \neq \emptyset. \qquad (18)$$

Being the PC expansion of the mapping $\xi \mapsto y(\xi)$ available i.e., the PC coefficient have been evaluated, the Sobol
indexes can be approximated as

$$S_u \approx \frac{\sum_{k \in K_u} c_k^2 \langle \Psi_k, \Psi_k \rangle}{\sum_{k=0}^{P} c_k^2 \langle \Psi_k, \Psi_k \rangle}, \quad \text{with} \quad K_u = \{k \in \{1, \dots, P\} | \Psi_k(\xi) = \prod_{i=1}^{|u|} \phi_{\alpha_i^k}(\xi_{u_i}), \alpha_i^k > 0\}, \quad (19)$$

where the symbol $\phi_k$ denotes the one-dimensional polynomial composing the orthogonal basis.

Being the mathematical formulation of the PC expansion established, we are left with the problem of determining the vector of PC coefficients $c$. Different approaches are available. Classically, a distinction among intrusive and non-intrusive methods is recognized. In general, intrusive methods rely on a Galerkin projection of the model equations on the PC basis and seek for the weak solution w.r.t. the random input. In this perspective, such technique requires access to the innermost part of the model which, in case of CFD applications, means accessing the source code of a vast and complex computational solver. On the other hand, non-intrusive methods only require a more or less extended set of deterministic realizations from the model under investigation. In the non-intrusive setting, the simulation framework is therefore employed as a mere black box query tool. A finite set $\xi^{(i)}$ of $d$ inputs, which can be sampled by means of different strategies, is associated to the corresponding output $y$ through the deterministic evaluation of the full computational model. Once a sufficient amount of realizations is available, the expansion coefficients of the linear regression can be reasonably approximated.

The procedure described above works quite well in most of the applications involving a smooth regular target function $f(x)$. Unfortunately, in the context of investigating icing on aeronautic airfoils, the highly non-linear dependence of the ice thickness model from the uncertain input variables is particularly difficult to capture, especially in the close proximity of the ice accretion limits. In such regions, the value of the predicted thickness lie on a constant lower bound (the zero value) for a vast portion of the input parameters space. In this perspective, the classical linear regression method described so far does not suit the approximation of continuous functionals achieving a constant lower bound on some part of their support. In the following, we formulate a possible solution, presenting a descriptive exposition of the idea underlying the proposed approach. Anticipating what will be carefully explained next, we propose a non-linear regression approach relying on the linear approximation of the targeted function, the ice thickness, composed with a thresholding function. For this reason, we refer to this formulation as a thresholded PC approximation. In the result Section IV, we will provide a comparison of the performances achieved with the standard linear approach against the ones achieved using the proposed non-linear regression formulation.

## A. Non-linear regression approach for in-flight icing applications

We propose a non-linear regression approach relying on the linear approximation of the targeted function, the ice thickness, composed with a thresholding function. Note that, since the resulting formulation is non-linear, the techniques generally employed for building standard (linear) PC approximations are no longer suitable or efficient. Therefore, in
our implementation we take advantage of a coordinate descent strategy to solve the minimization least square problem, to ultimately evaluate of the expansion coefficients. To enhance the stability and the quality of the regression, we also take advantage of the elastic net regularization strategy. The technical details of our implementation are thoroughly described in Appendix V. Besides, it is worth stressing that the thresholding function introduces a non-linearity resulting into the disruption of the analytic properties of standard (linear) PC. In order to characterize the QoIs statistics i.e., mean, variance and Sobol indexes, we employ surrogates in junction with a Monte Carlo sampling strategy.

**Problem statement** Consider $f: X \to \mathbb{R}$, where $X \subseteq \mathbb{R}^d$, such that $f(\mathbf{x}) = f_{\text{inf}}$ for $\mathbf{x} \in X_{\text{inf}} \subset X$ and $f(\mathbf{x}) > f_{\text{inf}}$ for $\mathbf{v} \in X \setminus X_{\text{inf}}$. Note that, without any loss of generality, we are not including any explicit reference to the time variable dependency $t$. We further assume that $f$ is continuous over $X$, that $|X_{\text{inf}}| > 0$ and $|X_{>}| > 0$, and that the lower bound $f_{\text{inf}}$ is known. We restrict ourselves to the case $f_{\text{inf}} = 0$; the generalization to $f_{\text{inf}} \neq 0$ is straightforward considering the approximation of $g(\mathbf{x}) = f(\mathbf{x}) - f_{\text{inf}}$.

We seek to construct an approximation of $f$ from a set of data consisting in $N$ couples $(x_i, f_i)$. Let $\mathcal{I} = \{1, 2, \dots, N\}$ be the set of data indices; we partition $\mathcal{I}$ into two distinct subsets: the set $\mathcal{I}_{\text{inf}}$ of data indexes with values $f_i = 0$, and $\mathcal{I}_{>} = \{i: f_i > 0\}$. The approximation is sought as a polynomial expansion composed with a thresholding function. Specifically, recalling that inputs are parametrized such that $\xi \doteq \xi(\mathbf{x})$, we set

$$f(\xi) \approx H_{\text{inf}} \left( \sum_{k=0}^{P} c_k \Psi_k(\xi) \right), \qquad (20)$$

where $\{\Psi_1, \dots, \Psi_P\}$ is a polynomial basis, $c_{k=1,\dots,P}$ are the unknown expansion coefficients and the thresholding function $H_{\text{inf}}$ is selected as

$$H_{\text{inf}}(u) = \begin{cases} u & \text{for } u \ge 0, \\ 0 & \text{for } u < 0. \end{cases} \qquad (21)$$

For notational convenience, we shall write $\mathbf{f} = (f_1 \cdots f_N)^\top \in \mathbb{R}^N$ the vector of observations, $\mathbf{c} = (c_1 \cdots c_P)^\top \in \mathbb{R}^P$ the vector of expansion coefficients, $[Z] \in \mathbb{R}^{N \times P}$ the matrix of regressors with entries $[Z]_{ik} \doteq \Psi_k(x_i)$ and $\tilde{\mathbf{f}} = [Z]\mathbf{c}$ the evaluation of the polynomial expansion $\tilde{f}(\xi) \doteq \sum_k c_k \Psi_k(\xi)$ at the data points (approximation before thresholding).

The objective is then to determine the vector of expansion coefficients $\hat{\mathbf{c}}$ that minimizes the approximation error. The function $f(\xi)$ being unknown, except at the data points, the minimization of the approximation error is substituted with a regression problem consisting in the generic minimization problem

$$\min_{\mathbf{c} \in \mathbb{R}^P} L(\mathbf{f} - H_{\text{inf}}([Z]\mathbf{c})), \qquad (22)$$

where $L$ is prescribed loss function and $H_{\text{inf}}$ applies to vectors component-wise. A classical choice for the loss is the
squared residual function, leading to the non-linear least squares minimization problem:

$$ \min_{\mathbf{c} \in \mathbb{R}^K} \|\mathbf{f} - H_{\text{Inf}}([Z]\mathbf{c})\|^2 = \sum_{i \in \mathcal{I}} [f_i - H_{\text{Inf}}([Z]\mathbf{c})]^2. \quad (23) $$

The least squares problem may require regularization to ensure its well-posedness. At this stage, we shall delay the discussion of regularization (which is reported in Appendix V) and focus on the solution of (23) in Sec. III.B.

## B. Alternative LS problem

A major difficulty in solving the LS problem (23) comes from the non-linear character and lack of differentiability of $H_{\text{Inf}}$. To avoid these difficulties, we propose an alternative formulation which does not explicitly rely on the thresholding function $H_{\text{Inf}}$. We start by observing that the summation in (23) can be split according to the subset of data indices; doing so it comes

$$ \begin{aligned} \sum_{i \in \mathcal{I}} (f_i - H_{\text{Inf}}(\tilde{f}(x_i)))^2 &= \sum_{i \in \mathcal{I}_{\text{Inf}}} H_{\text{Inf}}(\tilde{f}_i)^2 + \sum_{i \in \mathcal{I}_>} (f_i - H_{\text{Inf}}(\tilde{f}_i))^2 \\ &= \frac{1}{4} \sum_{i \in \mathcal{I}_{\text{Inf}}} (|\tilde{f}_i| + \tilde{f}_i)^2 + \sum_{i \in \mathcal{I}_>} (f_i - H_{\text{Inf}}(\tilde{f}_i))^2 \\ &= \frac{1}{4} \sum_{i \in \mathcal{I}_{\text{Inf}}} (2|\tilde{f}_i|^2 + 2|\tilde{f}_i|\tilde{f}_i) + \sum_{i \in \mathcal{I}_>} (f_i - H_{\text{Inf}}(\tilde{f}(x)))^2 \\ &\approx \frac{1}{2} \sum_{i \in \mathcal{I}_{\text{Inf}}} |f_i - \tilde{f}_i|^2 + \sum_{i \in \mathcal{I}_>} |f_i - \tilde{f}_i|^2 + \frac{1}{2} \sum_{i \in \mathcal{I}_{\text{Inf}}} |\tilde{f}_i|\tilde{f}_i \\ &= \sum_{i \in \mathcal{I}} w_i |f_i - \tilde{f}_i|^2 + \frac{1}{2} \sum_{i \in \mathcal{I}_{\text{Inf}}} |\tilde{f}_i|\tilde{f}_i. \end{aligned} \quad (24) $$

To obtain the final expression, we have used successively the relation $H_{\text{Inf}}(u) = |u| + u$ and the substitution of $H_{\text{Inf}}(\tilde{f})$ by $\tilde{f}$ for the data in $X \setminus X_{\text{Inf}}$. In other words, we consider the thresholding only at locations where $f_i = 0$ in the regression problem. Finally, the weight factors associated to the data are given by

$$ w_i = \begin{cases} 1 & \text{for } i \in \mathcal{I}_>, \\ 1/2 & \text{for } i \in \mathcal{I}_{\text{Inf}}. \end{cases} \quad (25) $$

As a result of these manipulations, the non-linear least squares problem reduces to a standard weighted least squares problem with an additional contribution whose treatment requires a special procedure, the gradient descent strategy, which is thoroughly described in Appendix V.
## IV. Results

The tests considered hereinafter aim at replicating experiments carried out at NASA's Glenn Icing Research Tunnel (IRT) facility [19, 20]. Each test involves a two-dimensional subsonic airfoil exposed to icing conditions within the cryogenic wind tunnel. The analyses reported in the following concern the propagation of the input aleatory uncertainty to the computational model outputs, by means of the UQ techniques introduced in Sec. III. PC surrogates of the computational model map the input parameters to the targeted QoI at discrete locations over the airfoil surface. Namely, a PC expansion is constructed for each element of the surface mesh, to characterize the numerical uncertainty at diverse locations along the airfoil curvilinear abscissa. Once the output statistics are characterized, numerical predictions, complemented by uncertainty estimates, are compared against experimental observation and repeatability data. The goal is to provide a preliminary accuracy assessment of the ice accretion prediction framework. Besides, performing a sensitivity analysis w.r.t. the input parameters helps shedding light on the icing mechanisms. For all the considered tests, a preliminary analysis (not reported) was carried out to define the proper numerical setting of computational simulations e.g., mesh resolution.

In the following, three different test cases are presented. In Section IV, we characterize the experimental uncertainty. The first test case, targeting the collection efficiency only, is presented in Sec. IV.A and concerns a GLC-305 airfoil. In the second test case, Sec. IV.B, a rime accretion experiment on a NACA0012 is replicated and the input uncertainty is propagated to the ice thickness layer predictions. Eventually, in Sec. IV.C we target a glaze ice accretion experiment over a NACA0012 profile.

### Characterization of the experiment

Ice accretion depends on a set of uncertain input parameters characterizing the environment and the flight conditions. Some of these are related to the aerodynamic field developing around the airfoil e.g., the Mach number (M), the angle of attack (AoA), the freestream pressure ($P_∞$) and temperature ($T_∞$), some other are related to the cloud e.g., the LWC and the MVD, whereas some others are related to the phase-change model e.g., the ice density $ρ_{ir}$ and the thermal conductivity $k_{ir}$.

We characterize input uncertainties by assuming that parameters are uniformly distributed within an interval centered on their reference nominal value. Case-by-case, the nominal value and the extent of the probability distribution support are imposed according to information reported in cited works. The only exception concerns the characterization of the uncertainty affecting the ice thermal conductivity $k_{ir}$, for which a slightly diverse approach is adopted. Indeed, $k_{ir}$ suffers from a dependency on the porous ice density (the air/ice volume ratio) and temperature. Through Eq. (11), the dependency of $k_{ir}$ on the rime ice density (and on some of the uncertain inputs) is made explicit. At the same time, variations of the porous ice thermal conductivity w.r.t. temperature correspond approximately to the one of pure ice [32] (at least in the range of air-to-ice ratios typical of in-flight icing). Because of this, the dependence of $k_{ir}$ on
temperature can be retrieved by applying a scaling factor, see Sec. II.C. The values of thermal conductivity of pure ice ($\hat{k}_i$) is available through tabulated values collecting measurements at different temperatures, see Tab. 1.

<table><caption>**Table 1** Pure glaze ice thermal conductivity tabulated w.r.t. temperature.</caption><tr><td>$T$ [K]</td><td>0.0</td><td>-5.0</td><td>-10.0</td><td>-15.0</td><td>-20.0</td><td>-25.0</td><td>-30.0</td><td>-40.0</td><td>-50.0</td><td>-55.0</td><td>-60.0</td></tr><tr><td>$\hat{k}_i$ [W/Km²]</td><td>2.22</td><td>2.25</td><td>2.30</td><td>2.34</td><td>2.39</td><td>2.45</td><td>2.50</td><td>2.57</td><td>2.63</td><td>2.76</td><td>2.90</td></tr></table>

To fit data and approximate $\hat{k}_i$ within the whole temperature range, a Gaussian process regressor with an exponentiated quadratic kernel function is employed [36]. A Karhunen-Loève decomposition of the posterior predictive covariance is applied, in order seek for a possible reduction of the dimensionality of the stochastic process. The spectral energy analysis of the posterior covariance reveals that the first mode contains most of the energy of the system. Therefore, the uncertainty affecting $\hat{k}_i$ is modeled through a single random variable $\xi_{\hat{k}_i} = N(0, 1)$, the Karhunen-Loève first mode. The reduced model is able to accurately fit the tabulated data up until a temperature of -60°C, which is well below the minimum value considered for the test cases presented hereinafter.

## A. Collection Efficiency

The first experiment concerns the GLC-305 airfoil exposed to a subsonic mist flow reproducing warm weather conditions. Temperature is well above the water freezing point and, therefore, the current investigation is focused on collection efficiency predictions only. The test-rig detailed set up and the complete measurements database are reported in [19]. Table 2 reports the experiment nominal conditions together with the uncertainty maximum and minimum bounds. Uncertainty bounds are established according to the IRT calibration report [37].

<table><caption>**Table 2** Test case A. Nominal test conditions and uncertainty bounds.</caption><tr><td>Mach</td><td>Pressure [Pa]</td><td>Temperature [K]</td><td>AoA [deg]</td><td>MVD [μm]</td></tr><tr><td>0.23 ± 0.0025</td><td>101325 ± 134</td><td>291.2 ± 0.5</td><td>1.6 ± 0.2</td><td>21 ± 1</td></tr></table>

The goal is to propagate the uncertainty affecting the parameters of the experiment namely, the angle of attach AoA, the mach number $M$, the free stream pressure and temperature $P_\infty$ and $T_\infty$ and the MVD, to the numerical predictions of the collection efficiency over the airfoil surface. For this test case, surrogates are constructed using the standard (linear) regression approach. A Legendre polynomial basis of order 3 is selected. The expansion coefficients are computed via the LS approach and 168 training points are employed. Additionally, 50 points where selected from the parameter space, using Latin Hyper-cube sampling, to define a test set of deterministic model realizations later exploited to assess the accuracy of surrogates.

Figure 5(a) reports the PCs mean collection efficiency (continuous line, left scale) and the root mean square error of the surrogates' prediction (dashed line, right scale) evaluated on the aforementioned test set. Figures are plotted
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Test case A. Collection efficiency ($\beta$) surrogate assessment. (a) Mean $\beta$ (left scale) and RMSE (right scale) of PC predictions evaluated on the test set. (b) Comparison among numerical predictions and experimental observations. Curves are complemented by the $2\sigma$ uncertainty obtained from PC of the computational model and experiment repeatability studies.</figcaption></figure>

against the curvilinear abscissa originating at the leading edge of the profile, having positive value in the direction of the pressure side. The RMSE is relatively small throughout the whole airfoil, thus proving the surrogates predictive capability for samples departing from the training ones.

Figure 5(b) reports the PCs mean collection efficiency (red curve), complemented by the $2\sigma$ uncertainty estimates (shaded red area), and the experimental observations (blue curve), complemented by the experiment repeatability envelope (shaded blue area). In addition, $\beta$ predictions from the NASA LEWICE solver from Ref. [19] are also reported for benchmark (black curve). The comparison reveals the predictive nature of our in-house solver, at least for the considered test conditions. The experiment repeatability region overlaps the computational solution variability almost everywhere. Differences, possibly due to the rough approximation of using the MVD parameter to characterize the cloud, are found near the impingement limits. A more accurate description of the cloud i.e., dropping the *MVD* parameter in favor of a full probability distribution of droplets diameters, could possibly help removing discrepancies w.r.t. the experiment.

In Fig. 6(a), we report the analysis of the $\beta$ variance ($\sigma_{\beta}^2$) for different PC truncation orders. Clearly, expansions truncated to the first, second and third order return comparable curves, with just small differences found in the close proximity of the impingement limits (the gray shaded area represents the region spanned by the impinging limits w.r.t. variation of the input parameters). On the left hand side, $\beta$ quickly grows from a zero value to a local maximum right past the end of the shaded area. This indicates that, despite the fact that the impinging limit remains confined within a quite limited region, the amount of water collected on the airfoil suction side suffers from a strong sensitiveness w.r.t. the input parameters variability. In other words, the region of the airfoil suction side in the close proximity of the stagnation point is critical to icing applications, in particular for what concern the design and the power requirements for
<figure><img src="image_7.png" /><figcaption>Fig. 6 Test case A. Collection efficiency ($\beta$) uncertainty analysis. (a) Collection efficiency variance for different PC truncation orders. The shaded area represents the region of the airfoil surface in to which the impingement limits are found with input parameters varying within the preset uncertainty bounds. (b) First order Sobol indexes.</figcaption></figure>

fully-evaporative anti-ice protection systems. On the right hand side i.e., on the pressure side, the curve shows a gentler variation leading to a quite large plateau which eventually decays to a local minimum, in the close proximity of the stagnation point, at $s \approx 0.002$. The variability of $\beta$ is limited to lower values but, interestingly, the region spanned by the suction side impinging limit is now wider. This points out another interesting fact which is (again) of the utmost relevance for the design of anti-ice fully-evaporative systems. Indeed, the portion of the surface to be protected is way larger than the one predicted from nominal condition. Therefore, designers should account for such a large displacement of the lower impinging limit to ensure safety even when flight conditions are off nominal.

Figure 6(b) shows instead the variance decomposition w.r.t the uncertain inputs i.e., first order Sobol indexes. According to this analysis, the AoA and the MVD uncertainties are found to be of the utmost relevance for icing problems since they provide the largest contribution to the $\beta$ variability, at least w.r.t. the considered (typical) uncertainty ranges. Three different zones on the airfoil can be highlighted. On the suction side, between the curvilinear coordinates $-0.017$ [m] and $-0.008$ [m] the variance shows a global maximum corresponding to one of the tails of the collection efficiency. Here, the variance is mostly due to uncertainties in the AoA ($\approx 80\%$), whereas the MVD provides a minor contribution ($\approx 20\%$). In correspondence with the peak of the collection efficiency, it is possible to identify a second zone in which the MVD uncertainty contributes the most (up to above 90%). Note that in this zone, the variance of the QoI reaches a local minimum. Continuing rightwards, we encounter a third region corresponding to the airfoil pressure side. Here, the AoA contributes by $\approx 80\%$ to the total variance whereas the remaining $\approx 20\%$ is again ascribable to the MVD (similarly to the first region, although the variance magnitude is here significantly lower). Since the associated Sobol indexes are negligible, there is a negligible contribution to the catch efficiency variance from the remaining uncertain parameters, Mach, pressure and temperature. To provide a means for quantifying the overall effect of each input, Tab. 3 reports the
<table><caption>Table 3 Test case A. Integral values, w.r.t. the curvilinear abscissa, of Sobol indexes normalized by the collection efficiency $\beta$ total variance</caption><thead><tr><th>Mach</th><th>Pressure</th><th>Temperature</th><th>AoA</th><th>MVD</th><th>higher order</th></tr></thead><tbody><tr><td>0.185%</td><td>0%</td><td>0.001%</td><td>28.92%</td><td>70.16%</td><td>0.72%</td></tr></tbody></table>

<table><caption>Table 4 Test case B. Nominal test conditions and uncertainty bounds.</caption><thead><tr><th>Mach</th><th>Pressure [Pa]</th><th>Temperature [K]</th><th>MVD [µm]</th><th>AoA [deg]</th><th>LWC [g/m³]</th><th>ξ<sub>k<sub>i</sub></sub></th></tr></thead><tbody><tr><td>0.21 ± 0.0025</td><td>100000 ± 134</td><td>244.8 ± 0.5</td><td>20 ± 1</td><td>4 ± 0.3</td><td>1 ± 0.07</td><td>N (0, 1)</td></tr></tbody></table>

integral value of each Sobol index w.r.t. the curvilinear coordinate. Results are normalized by the integral value of the $\beta$ total variance.

## B. Rime ice accretion

The aim of this test case is that of propagating the input uncertainties to ice thickness predictions. Computations reproduce a rime-ice accretion experiment performed at the IRT [20] and focuses on investigating the ice forming at the leading edge of a NACA0012 airfoil. The reference experiment includes a repeatability analysis targeting different values of the ice accretion parameters, thus providing uncertainty bounds complementing observations. The test duration is 360 s, conditions are reported in Tab. 4 together with their uncertainty characterization. We recall here that numerical simulations are obtained using the so-called one-shot setting.

Rime conditions were selected to assess the variability of the ice shape predictions resulting from the simpler accretion mechanism. Indeed, the supercooled droplets freeze immediately upon impact on the aircraft surface, basically implying the ice growth rate almost directly proportional to the amount of water collected (in other words, proportional to $\beta$), see (7). Nevertheless, the proportionality constant i.e., the inverse of the rime ice density ($1/\rho_{ir}$), is itself uncertain since the structure of the frozen layer will include air bubbles that remain trapped. As shown in Sec. II, it is possible to relate $\rho_{ir}$ to the uncertain input parameters via the Jones' model (10). In summary, compared to the test case A, we now consider the additional uncertainty on the LWC and $\xi_{k_i}$. Therefore, the dimensionality of the random inputs space for this test case is 6. A total number of samples (360, three times higher than the number of PC coefficient) is selected. Again, the linear regression approach is employed to build the PC surrogates of the full computational model.

Figure 7(a) reports the mean ice thickness distribution (continuous curve, left scale) together with the RMSE computed by comparing PC prediction against a testing data set (dashed curve, right scale). The mean thickness develops as a regular smooth curve along the airfoil curvilinear abscissa, not presenting any horn, as expected in a rime-type accretion problems. The RMSE is relatively low throughout the whole domain, but some loss of accuracy is seen near the accretion limits. This is related to the fact that the icing model consists in a continuous functional bounded by a lower constant value (the zero ice thickness) for which a polynomial approximation may not be suited. In the
<figure><img src="image_8.png" /><figcaption>Fig. 7 Test case B. (a) Mean ice thickness and RMSE for the linear and thresholded PCE surrogates for the rime case. (b) Comparison between experimental data and computed ice shape for the NACA0012 rime case.</figcaption></figure>

close proximity of the ice accretion limits, a large number of training data points (the vast majority of them, in fact) is associated with a zero ice thickness. That is, a large portion of the input parameter space is mapped to a plateau whereas the remaining is mapped to positive non-zero values. In view of this facts, the high error at the accretion limits is mainly due to linear PC surrogates failing in predicting a strictly non-negative value of the ice layer thickness. To overcome this issue, we took advantage of the non-linear formulation of the regression approach we propose in Sec. III. The RMSE of the PC thresholded approximation is reported in Fig. 7(a) for comparison w.r.t. the standard linear approach (dotted curve, right scale). The comparison reveals that the thresholded formulation generally enhances the accuracy of PC surrogates, in particular in the region of the ice accretion limits on airfoil pressure side (s ≈ 0.06).

Figure 7(b) compares the experimental ice shape reported in [20] with prediction from the PoliMIce framework. The area between the computed mean (red profile) plus/minus one standard deviation (red shaded area), is superimposed to the experimental envelope. The comparison of predictions and observations indicates the good predictive character of the simulation framework.

Figure 8(a) reports the ice thickness variance w.r.t. the location over the airfoil surface. Everywhere but near the bulk of the ice thickness (s ≈ 0.005), the variance approximately resembles that of the $σ_β^2$ curve. Similarities are due to the fact that, in the rime regime, impinging droplets freeze immediately upon impact. Therefore, the ice thickness mainly depends on the collection efficiency. Nevertheless, an additional variability peak is now present in correspondence of the stagnation point. Possibly, the presence of this additional peak is due to the rime ice thickness being linearly dependent on the LWC. This latter claim is supported by the analysis of variance, achieved by means of the Sobol decomposition, reported in Fig. 8(b). Clearly, the MVD, the AoA and the LWC uncertain parameters have a predominant role w.r.t. the variability of the targeted QoI. The AoA and the MVD dominate near the accretion limits whereas LWC is indeed of the utmost relevance in the region corresponding to the stagnation point. Instead, the contribution of the remaining uncertain
<figure><img src="image_9.png" /><figcaption>Fig. 8 Test case B. (a) Ice thickness variance w.r.t. the curvilinear abscissa. The shaded area represents the region of the airfoil surface in to which the ice accretion limits are found with input parameters varying within the preset uncertainty bounds. (b) First order Sobol indexes.</figcaption></figure>

<table><caption>Table 5 Test case B. Integral values, w.r.t. the curvilinear abscissa, of Sobol indexes normalized by the ice thickness (B) total variance.</caption><thead><tr><th>Mach</th><th>Pressure</th><th>Temperature</th><th>MVD</th><th>AoA</th><th>LWC</th><th>k<sub>ir</sub></th><th>higher order</th></tr></thead><tbody><tr><td>1.83%</td><td>0%</td><td>0.07%</td><td>22.44%</td><td>16.68%</td><td>57.3%</td><td>0%</td><td>1.68%</td></tr></tbody></table>

inputs is negligible. In particular, the uncertainty concerning the ice thermal conductivity is negligible throughout the whole domain. This was somehow expected since here we are targeting a rime type accretion which is governed by Eq. (7). Table 5 reports the relative contribution of each uncertain input to the variance of the QoI. As for the previous test case, the tabulated values are computed by integrating each Sobol indexes along the curvilinear coordinate. After normalizing by the integral along the airfoil of the QoI variance, the results in the table are shown as percentages.

### C. Glaze ice accretion

The third test case aims at propagating the input uncertainties to the ice thickness on a NACA 0012 airfoil, for a glaze type accretion regime. The overall exposition is 360 s. Given the highly non-linear nature of the glaze ice accretion phenomenon, the proposed thresholded PC approach is employed and compared against the standard (linear) PC approach. Analogously to Sec. IV.B, we rely on a least squares approach using a LHS strategy to independently define the PC training and test sets, which again include, respectively, 360 and 50 samples. The characterization of the stochastic, uniformly distributed, inputs is reported in Tab. 6.

<table><caption>Table 6 Test case C. Nominal test conditions and uncertainty bounds.</caption><thead><tr><th>Mach [-]</th><th>Pressure [Pa]</th><th>Temperature [K]</th><th>MVD [µm]</th><th>AoA [deg]</th><th>LWC [g/m³]</th><th>ξ<sub>k<sub>i</sub></sub></th></tr></thead><tbody><tr><td>0.23 ± 0.0025</td><td>100000 ± 134</td><td>267 ± 0.5</td><td>20 ± 1</td><td>4 ± 0.23</td><td>1 ± 0.07</td><td>N (0, 1)</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9 Test case C. (a) Mean ice thickness and RMSE for the linear and thresholded PCE surrogates. (b) Single random PC realization. Comparison of a single PC realization, obtained for the same input, from the linear and the thresholded PC surrogates.</figcaption></figure>

Figure 9(a) reports the mean ice thickness along the airfoil as a black continuous curve. The linear and the thresholded formulation result into overlapping results. The same plot reports, on the right axis, the RMSE for both approaches. By using the thresholded formulation, the error is again substantially reduced. Differently from test case B, here the most relevant enhancement is found in correspondence of the suction side accretion limit. The benefits of the thresholded formulation, w.r.t the linear one, can be better understood by comparing single realizations. Figure 9(b) reports the linear and the thresholded PC realizations of a point selected from the test set. At the same time, the reference solution obtained from the evaluation of the full computational model is also reported. The superior performance of the thresholded formulation can be clearly seen near the suction side accretion limit (s ≈ -0.015). At such locations, linear PC fails in predicting a zero ice thickness, returning negative, and thus not physical, values. On the contrary, the thresholded formulation is able to bound predictions to the zero thickness lower limit, recovering the accretion phenomenon with much grater accuracy and thus indicating superior performances.

Figure 10(a) reports the comparison of numerical prediction against experimental observations from [20]. The comparison reveals a fair agreement between simulations and measurements. Compared to the rime type accretion analyzed in test case B, we are now dealing with a much more irregular ice shape. The glaze accretion regime leads to the formation of a typical double-horned configuration. Moreover, the surface of the ice layer is now characterized by more pronounced ridges. The one-sigma uncertainty envelope associated to numerical predictions suffers from a quite sudden enlargement in the proximity of the accretion limits, suggesting how the test conditions uncertainty challenges the prediction of ice formation in such regions. To highlight this latter aspect, Fig. 10(b) reports the variance of numerical predictions due to the considered uncertain parameters, w.r.t. the airfoil curvilinear abscissa. A maximum can be found on the airfoil suction side, within the gray shaded region indicating the upper accretion limit range.

Figure 11(a) and (b) report, respectively, the first order Sobol indexes for $β$ and for $B$. By focusing on the ice
<figure><img src="image_11.png" /><figcaption>Fig. 10 Test case C. (a) Ice thickness distribution on NACA0012 airfoil. Computed and experimental bounds. (b) Ice thickness variance w.r.t. the curvilinear abscissa. The shaded area represents the region of the airfoil surface in to which the ice accretion limits are found with input parameters varying within the preset uncertainty bounds.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Test case C. First order Sobol indexes. (a) Collection efficiency $\beta$. (b) Ice layer thickness $B$.</figcaption></figure>
<table><caption>Table 7 Test case C. Relative contribution to the ice thickness variance throughout the whole airfoil surface.</caption><thead><tr><th>Mach</th><th>Pressure</th><th>Temperature</th><th>MVD</th><th>AoA</th><th>LWC</th><th>k<sub>ir</sub></th><th>higher order</th></tr></thead><tbody><tr><td>1.53%</td><td>0.02%</td><td>33.80%</td><td>10.96%</td><td>7.91%</td><td>23.73%</td><td>0.02%</td><td>22.04%</td></tr></tbody></table>

<table><thead><tr><th></th><th>μ [mm]</th><th>σ<sup>2</sup> [m<sup>2</sup>]</th></tr></thead><tbody><tr><td>s<sub>upper</sub></td><td>62.92</td><td>7.3948e-6</td></tr><tr><td>s<sub>lower</sub></td><td>-13.75</td><td>2.8905e-6</td></tr></tbody></table>

<table><caption>Table 8 Test case C. Mean, variance and normalized Sobol indexes for the accretion limits in the NACA0012 test case.</caption><thead><tr><th></th><th>σ<sup>2</sup><sub>Mach</sub> [%]</th><th>σ<sup>2</sup><sub>Pres</sub> [%]</th><th>σ<sup>2</sup><sub>Temp</sub> [%]</th><th>σ<sup>2</sup><sub>MVD</sub> [%]</th><th>σ<sup>2</sup><sub>AoA</sub> [%]</th><th>σ<sup>2</sup><sub>LWC</sub> [%]</th><th>σ<sup>2</sup><sub>K</sub> [%]</th><th>higher order [%]</th></tr></thead><tbody><tr><td>s<sub>upper</sub></td><td>0.4644</td><td>0.0133</td><td>0.0022</td><td>68.8246</td><td>30.3417</td><td>0.002998</td><td>0.01075</td><td>0.3399</td></tr><tr><td>s<sub>lower</sub></td><td>1.0863</td><td>0.0334</td><td>42.9033</td><td>10.3319</td><td>14.2846</td><td>30.3769</td><td>0.0362</td><td>0.9474</td></tr></tbody></table>

thickness, Fig. 11(b) shows that the uncertainty related to the free-stream temperature accounts for up to 90% of the total variance in the close proximity of the stagnation point. Indeed, in the glaze accretion regime, the freezing of the liquid film over the airfoil is strongly affected by the temperature profile. Near the accretion limit, where $\sigma_B^2$ is larger, the uncertainties on the LWC, the MVD and the AoA provide instead higher contributions. These parameters all concur to the determination of the trajectory of the particles (AoA and MVD) and of the incoming mass (LWC). Table 7 reports the overall contribution of each stochastic parameter to the total variance of B integrated over the whole airfoil. As seen earlier, the reported values are obtained by integrating the pointwise Sobol indexes along the curvilinear coordinate.

Eventually, since the prediction of the ice accretion limits is critical, it is interesting to further extend our UQ analysis to investigate this aspect. Therefore, the final step of our work targets the upper (and lower) icing limits position in terms of curvilinear abscissa coordinate. Table 8 reports the mean, variance and Sobol indexes of such QoIs.

For this particular test case, the position of the suction side accretion limit ($S_{upper}$) shows a variability larger than the pressure side one. Also, the relative contribution of the inputs to the total variance shows a great asymmetry. The suction side limit is mostly influenced by the MVD and AoA uncertainties which make up the 98% of the total variance, whereas the pressure side accretion limit is more influenced by the ice and water film dynamics. Indeed, on the lower side of the airfoil the variance of the icing limit is mainly influenced by the uncertainties in the free-stream temperature and LWC, with minor contributions from MVD and AoA.

<table><caption>Table 9 Test case C. Mean, variance and normalized Sobol indexes for the accretion limits in the NACA0012 test case. HO indicates higher order.</caption><thead><tr><th></th><th>μ [mm]</th><th>σ<sup>2</sup> [mm<sup>2</sup>]</th><th>σ<sup>2</sup><sub>Mach</sub></th><th>σ<sup>2</sup><sub>Pres</sub></th><th>σ<sup>2</sup><sub>Temp</sub></th><th>σ<sup>2</sup><sub>MVD</sub></th><th>σ<sup>2</sup><sub>AoA</sub></th><th>σ<sup>2</sup><sub>LWC</sub></th><th>σ<sup>2</sup><sub>K</sub></th><th>HO</th></tr></thead><tbody><tr><td>s<sub>upper</sub></td><td>62.92</td><td>7.395</td><td>0.4644%</td><td>0.0133%</td><td>0.0022%</td><td>68.82%</td><td>30.34%</td><td>0.002998%</td><td>0.01075%</td><td>0.3399%</td></tr><tr><td>s<sub>lower</sub></td><td>-13.75</td><td>2.891</td><td>1.086%</td><td>0.0334%</td><td>42.90%</td><td>10.33%</td><td>14.28%</td><td>30.37%</td><td>0.0362%</td><td>0.9474%</td></tr></tbody></table>
V. Conclusion and Future Work

A framework for simulating in-flight ice accretion under uncertainty is presented. For the first time, input uncertainties are propagated directly to the ice thickness through a full high-fidelity computational model, thus providing uncertainty bounds to the computed ice shapes, using a PC method. For all the presented cases, simulations fairly reproduce measurements collected at the NASA IRT. A comparisons of numerical predictions against experimental observations indicate the predictive capabilities of the computational framework.

Results expose and identify a major critical aspect consisting in a inherent difficulty in approximating the full icing model by means of PC surrogates. This aspect is particularly relevant in the close proximity of the ice accretion limits or, with a more general statement, at locations suffering ice build up only under particular, limited, circumstances. To mitigate such issue, a novel non-linear regression method for continuous functionals achieving a constant lower bound on some part of their support is proposed and exploited. Compared to the standard linear method, the proposed approach allows for the construction of PC surrogates endowed with superior accuracy performance.

Final results help identifying major sources of uncertainty and provide guidelines for the analysis of data. Concerning the prediction of the collection efficiency, the UQ analysis identified the dominant role of the cloud MVD and airfoil AoA uncertainty. An interesting outcome, of great relevance for the design of fully-evaporative anti-ice protection systems, is that the airfoil suction side presents a well defined region prone to water collection. Indeed, the associated impinging limit does not vary much w.r.t. the input uncertainty. At the same time, the amount of water collected suffers from a great variability. In other words, there is a well defined zone which collects a highly variable amount of water. This translates into specific design requirements namely, the extension of the fully-evaporative anti-ice system could be limited to a just well defined region but, at the same time, requiring a highly variable amount of power, depending on the flight conditions. On the contrary, on the airfoil pressure side we observe a quite opposite behavior. Indeed, the pressure side collects a less variable amount of water, and thus the anti-ice system can be fed with a more or less constant amount of power, but the critical region varies significantly with the flying conditions, thus requiring the protection of a larger area. With respect to the design of protection systems, similar conclusions can be drawn from the UQ analysis targeting the ice accretion. In the rime ice accretion regime, the variability of the predicted thickness is largely associated to the cloud MVD and airfoil AoA. This behavior, similar to the collection efficiency one, was expected since rime ice results from the immediate freezing of impinging droplets upon impact. Nevertheless, the rime ice thickness value shows an additional relevant sensitiveness w.r.t. the cloud LWC. Glaze ice accretion is instead abiding by slightly different rules which lead to slightly different results. Namely, the glace ice type accretion suffers from a strong sensitiveness w.r.t. the ambient temperature, particularly in the region of the stagnation point. Interestingly, the LWC and the MVD still retain a prominent role whereas the contribution of the uncertainty related to the AoA is basically reduced to negligible effects. This is a significant outcome for the purpose of designing anti-ice or de-ice on board systems. Indeed, results point out the necessity of considering a wide range of operating conditions, including both rime and glaze ice type, and that the
significance level of inputs uncertainty changes with the accretion regime.

Future work is expected to deliver a thorough Verification & Validation of the computational framework. Moreover, we plan to account for additional sources of uncertainty. Specifically, we intend to introduce a better detailed description of the cloud, relaxing the use of the MVD in place of a more representative probability distribution characterizing droplets' diameter. We are also looking towards accounting for the uncertainty inherent the reduced order models employed to estimate the droplet drag coefficient and walls interaction parameters. The above is expected to be relevant with particular reference to Supercooled Large Droplets (SLD) conditions encounters (Appendix O conditions). Besides, additional work needs to be done in order to assess the role of the uncertainty related to the aerodynamic flow field (e.g., random flow velocity fluctuations due to freestream turbulence). Concerning the ice accretion model, besides the characterization of the uncertain ice properties, questions still remain unanswered about the evaluation of the heat transfer contributions due to convection, kinetic release or evaporation. Naturally, the large number of uncertain factors influencing the in-flight ice accretion mechanism poses inherent barriers due to the well known curse of dimensionality. Therefore, an interesting research opportunity would be developing efficient and scalable methods for propagating uncertainty in in-flight ice applications.

## Appendix: non-linear regression for in-flight icing

In this Appendix, we report the technical details concerning the implementation of the coordinate descent strategy and of the regularization approach we employ to solve the least square problem (23).

**Coordinate descent strategy** For the minimization of

$$ \sum_{i \in I} w_i |f_i - \tilde{f}_i|^2 + \frac{1}{2} \sum_{i \in I_{\text{inf}}} |\tilde{f}_i - \tilde{f}_i|, \qquad (26) $$

we rely on a coordinate descent strategy. In this approach, the expansion coefficients $c_k$ are considered sequentially and optimized one at a time. From the current model estimation of the data, $\tilde{\mathbf{f}} = [Z]c$, we define for $k = 1$ to $M$ the partial model $\tilde{\mathbf{f}}^{\backslash k} = \tilde{\mathbf{f}} - [Z]e_k c_k$ ($e_k$ is the $k$-th canonical vector); $\tilde{\mathbf{f}}^{\backslash k}$ corresponds to the current estimation with $c_k$ set to zero. Letting $\mathbf{z}^k$ be the $k$-th column of $[Z]$, a more convenient expression for $\tilde{\mathbf{f}}^{\backslash k}$ is

$$ \tilde{\mathbf{f}}^{\backslash k} = \tilde{\mathbf{f}} - \mathbf{z}^k c_k. $$

Then, for fixed $k$, we seek for $\alpha \in \mathbb{R}$ such that $\tilde{\mathbf{f}}(\alpha) = \tilde{\mathbf{f}}^{\backslash k} + \mathbf{z}^k \alpha$ minimizes the loss:

$$ L(\tilde{\mathbf{f}}^{\backslash k}, \alpha) = \sum_{i \in I} w_i |f_i - (\tilde{f}_i^{\backslash k} + z_i^k \alpha)|^2 + \frac{1}{2} \sum_{i \in I_{\text{inf}}} |\tilde{f}_i^{\backslash k} + z_i^k \alpha|(\tilde{f}_i^{\backslash k} + z_i^k \alpha). \qquad (27) $$
To deal with the last summation, we proceed as follows. Assuming $z_i^k \neq 0$, we let $\alpha_i = -\tilde{f}_i^{\lambda^k}/z_i^k$; it comes for all $i \in \mathcal{I}_{\text{inf}}$

$$|\tilde{f}_i^{\lambda^k} + z_i^k \alpha|(\tilde{f}_i^{\lambda^k} + z_i^k \alpha) = |\tilde{f}_i^{\lambda^k} + z_i^k \alpha|^2 \text{sign}((\alpha - \alpha_i)z_i), \quad (28)$$

where

$$\text{sign}(u) = \begin{cases} 1 & \text{for } u \ge 0 \\ -1 & \text{for } u < 0. \end{cases} \qquad (29)$$

Finally, letting $\mathbf{r} = \mathbf{f} - \tilde{\mathbf{f}}^{\lambda^k}$, we obtain

$$L(\tilde{\mathbf{f}}^{\lambda^k}, \alpha) = \sum_{i \in \mathcal{I}} w_i |r_i - z_i^k \alpha|^2 + \frac{1}{2} \sum_{i \in \mathcal{I}_{\text{inf}}} |r_i - z_i^k \alpha|^2 \text{sign}((\alpha - \alpha_i)z_i^k) = \sum_{i \in \mathcal{I}} w_i^*(\alpha) |r_i - z_i^k \alpha|^2, \quad (30)$$

where

$$w_i^*(\alpha) = \begin{cases} (1 + \text{sign}(\tilde{f}_i^{\lambda^k} + \alpha z_i^k))/2 & \text{for } i \in \mathcal{I}_{\text{inf}}, \\ 1 & \text{for } i \in \mathcal{I}_{>}. \end{cases} \qquad (31)$$

The value $\alpha$ minimizing $L(\tilde{\mathbf{f}}^{\lambda^k}, \alpha)$ satisfies

$$\alpha = \frac{(\mathbf{z}^k)^{\top}[W^*(\alpha)]\mathbf{r}}{(\mathbf{z}^k)^{\top}[W^*(\alpha)]\mathbf{z}^k}, \qquad (32)$$

where $[W^*(\alpha)] = \text{diag}(w_1^*(\alpha), \dots, w_M^*(\alpha))$. The last expression suggests to compute $\alpha$ as a fixed point, through the iterative sequence:

$$\alpha^{t+1} = \frac{(\mathbf{z}^k)^{\top}[W^*(\alpha^t)]\mathbf{r}}{(\mathbf{z}^k)^{\top}[W^*(\alpha^t)]\mathbf{z}^k}. \qquad (33)$$

Finally, the updated value of $c_k = \lim_{t \to \infty} \alpha^t$. We observe that during the iterations, the vector of residuals remains constant $(\mathbf{r} = \mathbf{f} - \tilde{\mathbf{f}}^{\lambda^k})$ and only the diagonal matrix $[W^*(\alpha)]$ must be updated for $i \in \mathcal{I}_{\text{inf}}$ using (31). Algorithm 1 outlines the procedure for solving the non-linear regression problem.

**Regularized problem** For the stability and quality of the regression solution, it may be desirable to introduce a regularization when considering the minimization of (26). The regularization can becomes particularly crucial if $K$ becomes close to the cardinality of $\mathcal{I}_{>}$, that is when number of non-zero data becomes comparable or, even worst, less than the number of coefficients to estimate. Typically, $\ell_1$ (Lasso), $\ell_2$ (Tychonoff), or a combination of both (elastic net)
**Algorithm 1** Outline of the iterative procedure (coordinate descent) to compute the expansion coefficients.

**procedure** FRT(Obsrvations **f**, matrix of regressor [Z], Initial coefficients **c**)

for all $i \in I$ do

Set $\mathbf{r} = \mathbf{f} - [Z]\mathbf{c}$

▷ initialize residual

end for

while Not converged do

for $k = 1, \dots, K$ do

Set $\mathbf{z} = [Z]_k$

$\mathbf{r} \leftarrow \mathbf{r} + \mathbf{z}c_k$

$\alpha = 0$

Set $[W^*](\alpha)$

while Not converged do

$\alpha \leftarrow \frac{\mathbf{z}^\top[W^*(\alpha)]\mathbf{r}}{\mathbf{z}^\top[W^*(\alpha)]\mathbf{z}}$

update $[W^*](\alpha)$

end while

Set $c_k = \alpha$

Set $\mathbf{r} \leftarrow \mathbf{r} - \mathbf{z}c_k$

end for

end while

end procedure

▷ Loop over the coordinates

▷ Set regressor to k-th column of [Z]

▷ Set residual for $c_k = 0$

▷ Initialize correction

▷ See (31)

▷ Iteration on correction

▷ Update correction

▷ Update the weights with (31)

▷ Update coefficient

▷ Update residual

are employed, extending (26) as follows:

$$ \sum_{i \in I} w_i |f_i - \tilde{f}_i|^2 + \frac{1}{2} \sum_{i \in I_{\text{inf}}} |\tilde{f}_i| |f_i + \beta \|\mathbf{c}\|_{\ell_2}^2 + \lambda \|\mathbf{c}\|_{\ell_1}, \quad (34) $$

with regularization parameters $\beta, \lambda \ge 0$.

The derivation presented above easily accommodate the regularization terms. The procedure remains the same, except that the optimum $\alpha$ now minimizes

$$ L_{\text{Reg}}(\tilde{\mathbf{f}}^k, \alpha) = \sum_{i \in I} w_i^*(\alpha) |r_i - z_i^k \alpha|^2 + \beta \alpha^2 + \lambda |\alpha|. \quad (35) $$

The solution, known as the soft-thresholded estimate, is given by

$$ \alpha = \text{sign}(\alpha_{\ell_2}) \max(0, |\alpha_{\ell_2}| - \lambda), \quad \alpha_{\ell_2} = \frac{(\mathbf{z}^k)^\top [W^*](\alpha) \mathbf{r}}{(\mathbf{z}^k)^\top [W^*](\alpha) \mathbf{z} + \beta}. \quad (36) $$

In the previous equation, we have denoted $\alpha_{\ell_2}$ the classical estimate for the Tychonoff regularization (note the $\beta$ terms that shrinks the coefficients). The coordinate descent algorithm 1 can be easily extended to account for the regularization.

As for the determination of the regularization parameters, classical procedures involve cross-validation. For instance, one can determine the regularization constant $\beta$ and $\lambda$ that minimize the mean squared prediction error (of the
thresholded approximation) with respect to an auxiliary validation set with data $\mathbf{f}'$ and regression matrix $[Z']$, that is

$$CV(\lambda, \beta) = \|\mathbf{f}' - H_{\text{inf}}([Z']\mathbf{c}(\beta, \lambda))\|_{2}^{2}. \qquad (37)$$

## Funding Sources

G. Gori would like to acknowledge that this work was partially funded by the European Commission's H2020 program through the UTOPIAE Marie Curie Innovative Training Network (H2020-MSCA-ITN-2016, Grant Agreement number 722734). A. Guardone acknowledges funding from the European Union's Horizon 2020 under grant agreement No. 824310 ICE GENESIS.

## References

[1] Stefan, J., "Ueber die theorie der eisbil-dung im polarmeere," *Annalen Der Physik Und Chemie*, Vol. 42, 1891, pp. 269-286.

[2] Messinger, B. L., "Equilibrium Temperature of an Unheated Icing Surface as a Function of Air Speed," *Journal of the Aeronautical Sciences*, Vol. 20, No. 1, 1953, pp. 29-42. https://doi.org/10.2514/8.2520, URL https://doi.org/10.2514/8.2520.

[3] Aliaga, C. N., Aubé, M. S., Baruzzi, G. S., and Habashi, W. G., "FENSAP-ICE-Unsteady: Unified In-Flight Icing Simulation Methodology for Aircraft, Rotorcraft, and Jet Engines," *Journal of Aircraft*, Vol. 48, No. 1, 2011, pp. 119-126.

[4] Beaugendre, H., Morency, F., and Habashi, W. G., "FENSAP-ICE's Three-Dimensional In-Flight Ice Accretion Module: ICE3D," *Journal of Aircraft*, Vol. 40, No. 2, 2003, pp. 239-247.

[5] Beaugendre, H., Morency, F., and Habashi, W. G., "FENSAP-ICE's Three-Dimensional In-Flight Ice Accretion Module: ICE3D," *Journal of Aircraft*, Vol. 40, No. 2, 2003, pp. 239-247. https://doi.org/10.2514/2.3113, URL https://doi.org/10.2514/2.3113.

[6] Gori, G., Zocca, M., Garabelli, M., Guardone, A., and Quaranta, G., "PoliMIce: A simulation framework for three-dimensional ice accretion," *Applied Mathematics and Computation*, Vol. 267, 2015, pp. 96-107. https://doi.org/10.1016/j.amc.2015.05.081, URL http://www.sciencedirect.com/science/article/pii/S0096300315007055, the Fourth European Seminar on Computing (ESCO 2014).

[7] Gori, G., Parma, G., Zocca, M., and Guardone, A., "Local Solution to the Unsteady Stefan Problem for In-Flight Ice Accretion Modeling," *Journal of Aircraft*, Vol. 55, No. 1, 2018, pp. 251-262. https://doi.org/10.2514/1.C034412, URL https://doi.org/10.2514/1.C034412.

[8] Hedde, T., and Guffond, D., "ONERA Three-Dimensional Icing Model," *AIAA Journal*, Vol. 33, No. 6, 1995, pp. 1038-1045.

[9] Mingione, G., and Brandi, V., "Ice Accretion Prediction on Multielement Airfoils," *Journal of Aircraft*, Vol. 35, No. 2, 1998, pp. 240-246.
[10] Myers, T. G., and Charpin, J. P. F., “A mathematical model for atmospheric ice accretion and water low on a cold surface,” *International Journal of Heat and Mass Transfer*, Vol. 47, 2004, pp. 5483–5500.

[11] Myers, T. G., “Extension to the Messinger Model for Aircraft Icing,” *AIAA Journal*, Vol. 39, No. 2, 2001, pp. 211–218.
https://doi.org/10.2514/2.1312, URL https://doi.org/10.2514/2.1312.

[12] Petrosino, F., Mingione, G., Carozza, A., Gilardoni, T., and D'Agostini, G., “Ice Accretion Model on Multi-Element Airfoil,” *Journal of Aircraft*, Vol. 48, No. 6, 2011, pp. 1913–1920.

[13] Verdin, P., Charpin, J. P. F., and Thompson, C. P., “Multistep Results in ICECREMO2,” *Journal of Aircraft*, Vol. 46, No. 5, 2009, pp. 1607–1613.

[14] W. B. Wright, D. G., R. W. Gent, “DRA/NASA/ONERA Collaboration on Icing Research Part II–Prediction of Airfoil Ice Accretion,” NASA Contractor Report 202349, NASA, 1997.

[15] Wright, W., “User’s Manual for LEWICE Version 3.2,” Contractor Report NASA/CR-2008-214255, NASA, 2008.

[16] Gent, G., Dart, N., and Cansdale, J., “Aircraft icing,” *Phil. Trans R Soc Lond A 358:2873 2911*, 2000.

[17] A. M. DeGennaro, L. M., C. W. Rowley, “Uncertainty Quantification for Airfoil Icing Using Polynomial Chaos Expansions,” *Journal of Aircraft*, Vol. 52, No. 5, 2015, pp. 1404–1411.

[18] L. Prince Raj, R. S. M., K. Yee, “Sensitivity of ice accretion and aerodynamic performance degradation to critical physical and modeling parameters affecting airfoil icing,” *Aerospace Science and Technology*, Vol. 98, 2020.

[19] Papadakis, M., Wong, S., Rachman, A., Hung, K., Vu, G., and Bidwell, C., “Large and small droplet impingement data on airfoils and two simulated ice shapes,” *NASA Technical Memorandum*, 2007, pp. 2007–213959.

[20] Shin, J., and Bond, T. H., “Results of an Icing test on a NACA 0012 airfoil in the NASA Lewis Icing Research Tunnel,” *NASA Technical Memorandum*, 1992, pp. 1992–105374.

[21] Le Maître, O., and Knio, O., *Spectral Methods for Uncertainty Quantification*, 1st ed., Scientific Computation, Springer Netherlands, 2010. https://doi.org/10.1007/978-90-481-3520-2.

[22] Palacios, F., Alonso, J., Duraisamy, K., Colonno, M., Hicken, J., Aranake, A., Campos, A., Copeland, S., Economon, T., Lonkar, A., Lukaczyk, T., and Taylor, T., *Stanford University Unstructured (SU&lt;sup&gt;2&lt;sup&gt;): An open-source integrated computational environment for multi-physics simulation and design*, 2013. https://doi.org/10.2514/6.2013-287, URL https://arc.aiaa.org/doi/abs/10.2514/6.2013-287.

[23] Menter, F., “Zonal Two Equation $k-\omega$ Turbulence Models for Aerodynamic Flows,” *AIAA Paper*, 1993.

[24] Crowe, C. T., Schwarzkopf, J. D., Sommerfeld, M., and Tsuji, Y., *Multiphase Flows with Droplets and Particles, 2nd ed.*, Engineering & Technology, Physical Sciences, Boca Raton: CRC Press, 2011. https://doi.org/10.1201/b11103.
[25] Loth, E., "Numerical Approaches for Motion of Dispersed Particles, Droplets and Bubbles," *Progress in Energy and Combustion Science*, Vol. 26, 2000.

[26] Morrison, A., *An Introduction to Fluid Mechanics*, Cambridge University Press, New York, 2013.

[27] R. Honsek, W. H., "Eulerian modeling of droplet impingement in the SLD regime of aircraft icing," *AIAA Paper*, 2006.

[28] R. Clift, M. W., J.R. Grace, *Bubbles, drops, and particles*, Courier Corporation, 2005.

[29] Bai, C., and Gosman, A., "Development of Methodology for Spray Impingement Simulation," Tech. rep., SAE Technical Report 950283, 1995.

[30] T. Bellosta, A. G., G. Parma, "A Robust 3D Particle Tracking Solver for in-Flight Ice Accretion Using Arbitrary Precision Arithmetic," *Proceedings of the VIII International Conference on Computational Methods for Coupled Problems in Science and Engineering*, edited by M. P. E. Onate and B. Schrefler, CIMNE, 2019, pp. 41-50.

[31] Rios, M. A., "Icing Simulations Using Jones' Density Formula for Accreted Ice and Lewice," *29th Aerospace Sciences Meeting*, 1991.

[32] N. Calonne, e. a., "Thermal Conductivity of Snow, Firn, and Porous Ice From 3-D Image-Based Computations," *Geophysical Research Letters*, Vol. 46, 2019.

[33] Reagan, M. T., Najm, H. N., Debusschere, B. J., Le Maître, O. P., Knio, O. M., and Ghanem, R. G., "Spectral stochastic uncertainty quantification in chemical systems," *Combustion Theory and Modelling*, Vol. 8, No. 3, 2004, pp. 607-632.
https://doi.org/10.1088/1364-7830/8/3/010, URL https://www.tandfonline.com/doi/abs/10.1088/1364-7830/8/3/010.

[34] Crestaux, T., Le Matre, O., and Martinez, J.-M., "Polynomial chaos expansion for sensitivity analysis," *Reliability Engineering & System Safety*, Vol. 94, No. 7, 2009, pp. 1161-1172.

[35] Sobol, I. M., "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates," *Mathematics and Computers in Simulation*, Vol. 55, No. 1, 2001, pp. 271-280.
https://doi.org/10.1016/S0378-4754(00)00270-6, URL http://www.sciencedirect.com/science/article/pii/S0378475400002706, the Second IMACS Seminar on Monte Carlo Methods.

[36] Rasmussen, C. E., and Williams, C. K. I., *Gaussian Processes for Machine Learning*, 2nd ed., MIT Press, 2006.

[37] E. A. Arrington, D. W. S., M. T. Pickett, "Flow Quality Studies of the NASA Lewis Research Center Icing Research Tunnel Diffuser," NASA Technical Memorandum 106311, NASA, 1994.