A rotorcraft in-flight ice detection framework using computational aeroacoustics and Bayesian neural networks

This work develops a novel ice detection framework specifically suitable for rotorcraft using computational aeroacoustics and Bayesian neural networks. In an offline phase of the work, the acoustic signature of glaze and rime ice shapes on an oscillating wing are computed. In addition, the aerodynamic performance indicators corresponding to the ice shapes are also monitored. These performance indicators include the lift, drag, and moment coefficients. A Bayesian neural network is subsequently trained using projected Stein variational gradient descent to create a mapping from the acoustic signature generated by the iced wings to predict their performance indicators along with quantified uncertainty that is highly important for time- and safety-critical decision-making scenarios. While the training is carried out fully offline, usage of the Bayesian neural network to make predictions can be conducted rapidly online allowing for an ice detection system that can be used in real time and in-flight.


Introduction
Rotorcraft possess the capability to vertically take-off and land, allowing them to operate in challenging flight conditions where conventional fixed-wing aircraft cannot.Their operational environment lends itself to low-level altitudes, and many missions are often toward the boundary of the flight envelope.Throughout the course of aviation, the requirement for flight in poor weather conditions and at night has always been a driver for innovation and rotorcraft are no exception.In the past, rotorcraft were fair weather vehicles with marginal performance; however, now they regularly operate in conditions from hot and dry to cold, wet, and windy (Padfield 2008).In particular, cold climates can potentially leave aircraft susceptible to dangerous icing conditions.In-flight icing encounters can jeopardize the performance and handling qualities of aircraft and hence pose a serious threat to flight safety (Bragg et al. 1986;Gent et al. 2000).This threat to flight safety arises as ice accretion can rapidly alter aerodynamic lifting surfaces such as wings and rotors during flight which are highly sensitive to geometric modifications.Naturally, establishing a comprehensive understanding of the fundamental physics of aircraft icing is paramount to ensure the reduction of ice-related accidents and the progression of safety-critical technology.
In 2002, Bragg et al. (2022) critically assessed aircraft safety during icing conditions and introduced smart icing systems to reduce aircraft accident rates.Their work recognized that warning the pilot of ice conditions through the use of ice detection systems could significantly increase the safety of aircraft to either optimize the use of ice protection systems or to avoid icing conditions altogether.Over the past two decades, this has led to sustained technological advancements in ice detection systems.
Methods for ice detection on aircraft can be characterized as either direct or indirect.The former involves using sensors to detect ice accretion on the airframe or atmospheric conditions symptomatic of icing environments, while the latter relies on monitoring for a change in response to the aircraft state due to icing.Direct methods have received the most scientific interest where many novel sensing techniques have been developed.One prominent direct strategy entails measuring the reaction of a signal or electric field to the presence of ice, such as through a surface-generated wave, infrared thermography, optical reflectometry, or ultrasonics wave techniques.For example, Varadan et al. (1997) proposed producing surface-generated Love waves from piezoelectric devices and propagating them through a surface.As ice accretes, the amplitude and velocity of these waves would be altered.Hongerholt et al. (2002) propagated and reflected ultrasonic-guided waves on the surface of a wing and measured the attenuated amplitude of the reflected waves in different mediums for detecting the existence of ice.Zhuge et al. (2012) used infrared thermography where emitted electromagnetic signals reflected by the surface are analyzed for temperature variations.Bassey and Simpson (2007) employed optical time domain reflectometry to measure the dielectric permittivity and volumetric water content on the aircraft surface.Schlegl et al. (2015) further used capacitive techniques to monitor the presence of dielectric substances, leveraging the different electrical properties of air, water, and ice.Roy et al. (1998) generated an excitation signal and measured the impedance change in the electrodes when ice is formed.
Overall, direct ice detection techniques allow for precise information regarding the ice location, thickness, and location.However, they are largely commercially unproven due to preliminary testing not exceeding the laboratory environment or are unpractical for moving components, making their implementation on rotorcraft challenging.
Indirect ice detection approaches are based on comparing the aircraft state to previously established clean states and using the difference between the iced and clean aircraft states as an indicator for icing.Melody et al. (2000) utilized aircraft dynamics as an indicator for icing, which required algorithms to identify changes in aircraft performance, stability, and control parameters over time based on measurements of the aircraft state variables and control input.Deiler and Fezans (2020) introduced a performance-based ice detection methodology based on the change in the steady flight state, which can be used prior to reaching dynamic stall conditions and during steady flight conditions contrary to the aircraft dynamics approach.However, a major limitation of using either the dynamics or performance of aircraft for ice detection is that the type, location, and thickness of ice cannot be determined.Additionally, the detection of ice prior to dynamic stall is fundamental to aircraft safety.
Another major type of indirect ice detection involves using acoustic noise.The acoustics-based techniques are particularly relevant for rotorcraft, which bear a characteristic and largely undesirable feature that is their rotor noise.Rotor noise sources are complex and are a combination of loading noise, thickness noise, shock noise, blade-vortex interaction noise, and broadband noise.However, Cheng et al. (2018) recognized the strong correlation between the ice-induced surface roughness and the broadband noise level suggesting rotor acoustics could be a viable technology for rotorcraft ice detection.Chen et al. (2019) further investigated using numerical techniques and an artificial iced notch and concurred that the acoustics could be potentially used to indicate the formation of ice on a rotor.The experimental work from Cheng et al. (2018) used sandpaper to represent ice roughness, while the numerical work from Chen et al. ( 2019) used an artificial ice notch based on the experimental iced surface roughness results of a straight wing from the University of Illinois, presented by Broeren and Bragg (2005).The first characterization between glaze and rime iced noise signals was established by Morelli et al. (2020).Their work utilized numerical techniques to simulate the accretion of ice on rotors and subsequently the rotor noise signature.The results highlighted that glaze ice accretion produces noise at a significantly higher frequency caused by ice-induced flow separation.In the main, ice detection via means of aircraft acoustics is still in the early stages of research.
This paper seeks to present a proof-of-concept acousticsbased ice detection system for rotorcraft through building a machine learning (ML) model.While several studies mentioned above have investigated using noise as an method for ice detection, this work is the first to develop a system that predicts the aerodynamic performance (e.g., lift and drag coefficients) directly from observed acoustics.In other words, we seek to not only state whether icing is present or not, but also predict the severity and consequence of the icing condition through the aerodynamic performance metrics directly useful for aviation and piloting.To achieve this, we first utilize numerical techniques to simulate a large dataset of ice shapes.Next, high-fidelity scale-resolving methods are used to solve for the flow field over the iced rotor blades and extract their corresponding aerodynamic performance indicators.The surface pressure fluctuations are then taken from the aerodynamic data and used by an aeroacoustic solver to predict the far-field noise levels.Finally, we train a ML model to directly map the acoustic noise signal produced from the ice structures to the aerodynamic performance indicators of the aircraft; this "shortcut" can allow Page 3 of 17 197 predictions to be made at real-time speeds.An overview of our overall framework is shown in Fig. 1.
We choose neural network (NN) as our ML model for its expressiveness in capturing complex non-linear mappings and in particular, Bayesian neural network (BNN) that also provides uncertainty quantification (UQ) in the ML model and its predictions.Other ML model types may also be adopted: for example, Gaussian process is another popular choice with built-in UQ, but it lacks scalability (covariance matrix inversions become expensive when more training data become available) and flexibility (typically limited for predicting low-dimensional quantities).In contrast, NNs are known for their expressiveness and ability to handle large datasets.Furthermore, the BNN framework can flexibly couple with advanced NN architectures such as convolutional neural networks (CNN) (Lecun et al. 1998) and long short-term memory units (LSTM) (Hochreiter and Schmidhuber 1997) that are specialized in high-dimensional data and time series signals, respectively.Indeed, NNs have already been used for in-flight ice detection and identification (Caliskan and Hajiyev 2013) albeit not with rotorcraft acoustics.We further highlight the critical need to provide UQ for the ML models used in such mission-and safety-critical settings-we need to know the quality of a ML prediction before using it for decision-making.However, performing UQ for a NN is a challenging task due to its high number of model parameters.A further contribution of our work is then to enable UQ computation in BNNs leveraging recently developed algorithm of projected Stein variational gradient descent (pSVGD) (Chen and Ghattas 2020).
Several limitations of this complex problem do however need to be acknowledged.While this methodology is able to predict the tonal and broadband noise components of iced airfoils, in reality the noise of an operating helicopter is highly complex with many different components including engine noise and tail rotor noise which would likely pollute the small differences in noise observed by an iced rotor.
Nevertheless, this is the first step as to understanding if an acoustic ice detection system using NN is feasible.To that end, if successful the installation on helicopter in operation would be possible due to the system being non-intrusive.
In summary, we develop a framework for rotorcraft icing detection highlighted by the following key novelty and contributions: • Uses in-flight acoustic measurements to enable prediction of aerodynamic performance for rotorcraft, who cannot use visual/image observations due to the high-blade RPM during operation; • Enables real-time performance predictions with quantified uncertainty through advanced ML models of BNNs, using a recently developed pSVGD algorithm; • Uses simulation data that are generated from state-of-theart physics-based models rather than empirical models; and • Requires a truly multi-disciplinary approach to combine ice formation and growth, turbulent flow aerodynamics, aeroacoustics, and ML with UQ.
The structure of the remaining paper is as follows: Sect. 2 introduces the methodology of the acoustic ice detection and the ML using BNNs; Sect. 3 discusses the results of the ice prediction, the turbulent flow simulations and acoustic prediction, as well as the training of the BNN and the design of the ice detection system; and Sect. 4 highlights the concluding remarks of the work.

Methodology
The development of a technique for real-time, in-flight aeroacoustic ice detection system requires a highly inter-disciplinary approach.The methodology used for the prediction of in-flight ice accretion, the simulation of the turbulent iced flow field and noise prediction, and finally the ML are hereinafter introduced.

Icing simulation
The PoliMIce ice prediction toolkit is utilized to simulate a large sample of ice shapes required for the ML model training.The PoliMIce ice prediction framework is described in detail in Gori et al. (2015a) and Morelli and Guardone (2021) and will subsequently be summarized.A conventional icing simulation structure routinely involves a three-stage process which iteratively updates to account for unsteady ice accretion.This simulation process is commonly known as multi-step ice accretion and is adopted by PoliMIce in this work.A schematic of multi-step ice accretion is shown in Fig. 2. The first stage of the icing framework is to simulate the aerodynamic flow field around the region of interest.In this work, the SU2 solver is used to compute the aerodynamic flow field over the wing, with further details of the SU2 solver discussed in Sect.2.2.
The second stage is to simulate the trajectories of supercooled water droplets within the flow domain and determine the collection efficiency and impingement locations.An in-house particle tracking code based on a Lagrangian framework was developed at Politecnico di Milano and is used for the simulation of clouds containing supercooled water droplets (Bellosta et al. 2019).The Lagrangian framework permits the modeling of supercooled water droplets at an individual particle level.Subsequently, fundamental physics such as particle splashing on impact can be modeled.Techniques to allow for particle tracking in mesh with arbitrary motion and across non-conformal boundary interfaces are used within this work to simulate clouds entrained within rotorcraft flow fields (Morelli et al. 2019).
Finally, an icing solver is utilized to predict the rate of ice accretion.The in-house code PoliMIce is used for computing the ice accretion (Gori et al. 2015a).The PoliMIce software library provides state-of-the-art ice formation models, including the Myers model (Myers 2001), the modified Myers model (Gori et al. 2015b), and the local exact solution of the Stefan problem (Gori et al. 2017).Myers introduced a rime limit thickness, B g , as a criterion for the selection of the ice regime.The rime ice limit thickness describes the condition at which the glaze regime can first appears.If the ice thickness B < B g or B g < 0 ∶ the rime ice accretion law is used, while if the ice thickness B > B g ∶ the glaze ice accretion law is used.The rate of rime ice accretion can then be written as while the rate of glaze ice accretion reads where the collection efficiency, droplet liquid water content, and freestream velocity are, respectively, denoted by , LWC , and U ∞ .Here, represents the density of ice which depends on the ice regime.Qdown and Qup are the heat fluxes exchanged on the ice/water phase change interface.L F is the latent heat required for melting ice.
The model used in this work to capture the complex experimental ice shapes is the local exact solution of the unsteady Stefan problem for the temperature profiles within the ice layer in glaze conditions (Gori et al. 2017).The ice shapes are then computed using a multi-step approach, whereby non-linear ice accretion is accounted for by iteratively updating the surface solution on which the ice accretes.The efficacy of the approach for predicting ice shapes on two-dimensional oscillating airfoils is highlighted in Morelli et al. (2020).

Turbulent flow simulation and noise prediction
SU2 has been developed with the task of solving partial differential equations (PDE) and constrained optimization problems on general unstructured meshes.The core of the suite is a Reynolds-averaged Navier-Stokes (RANS) solver capable of simulating the compressible, turbulent flows that are representative of many problems in aerospace and mechanical engineering.The finite volume method (FVM) is applied on arbitrary unstructured meshes using a standard edge-based data structure on a dual grid with control volumes constructed using a median-dual, vertex-based scheme.Several numerical fluxes like Jameson-Schmidt-Turkel (JST), Upwind Roe, and variants of the AUSM schemes are implemented and slope limiters enable second-order space integration.For time discretization SU2 uses a second-order dual-time stepping method with an outer time loop to march through the physical time and of an inner loop which is usually a pseudo-time iteration or a (quasi-)Newton scheme.An Arbitrary Lagrangian-Eulerian (ALE) formulation of the RANS equations is implemented in SU2 to account for unsteady grid motion, in which the convective fluxes are adjusted to obtain solutions on arbitrarily moving grids.
In the SU2 suite, of particular relevance to this work are the scale-resolving and acoustic prediction capabilities.For sub-grid scale model (Shur et al. 2015) was used.In addition, to limit the numerical dissipation in the LES part of the EDDES solver, the inviscid flux is computed using the Simple Low Dissipation Advection Upstream (SLAU2) scheme (Kitamura and Hashimoto 2016).The EDDES-SA solver has been demonstrated to successfully predict separated flows (Molina 2015;Molina et al. 2019), especially those induced by icing on the wing leading edge (Molina et al. 2020) For aeroacoustic prediction, a solid-surface Ffowcs Williams and Hawkings (FWH) analogy is used.To that end, Farassat's Formulation-1A (F1A) (Farassat 2007) has been implemented in the SU2 suite (Zhou et al. 2017;Icke et al. 2020;Morelli et al. 2022) which computes far-field noise from surface pressure fluctuations extracted from a preceding unsteady EDDES-SA simulation.

Machine learning using Bayesian neural networks
In this section, we first introduce the classical NN (Rosenblatt 1958), followed by the BNN (Bate et al. 1998) that captures the uncertainty of NN model parameters and finally the pSVGD (Chen and Ghattas 2020) algorithm used to numerically constructs a BNN.

Neural networks
A NN is a function f mapping input feature x to output target variable y.We write ŷ = f (x) where the hat denotes the predicted value from the NN.For illustration, we focus solely on densely connected feedforward NNs (also known as multi-layer perceptrons), although more advanced NN architectures (e.g., convolutional and recurrent NNs) may also be adopted within our overall ML framework.The hyperparameters that define the NN architecture (e.g., number of layers, nodes per layer, and activation functions) are selected prior to NN training, and good hyperparameter choices can be sought systematically through validation or cross-validation.Once the hyperparameters are chosen, the NN training then focuses on determining the weight and bias parameters of all NN layers.We denote the collection of all such trainable parameters by .Then, more explicitly, we write the prediction under a particular NN parameter setting as ŷ = f (x;) Given N training points in the form of input-output pairs (x T , y T ) = {x n , y n } N n=1 , NN training seeks * that minimizes a loss function reflecting the degree of mismatch between NN predictions and the true target values for these training points.In regression problems, for example, the mean squared error (MSE) loss is often used: This optimization problem is typically solved via gradient-based algorithms, such as stochastic gradient descent (LeCun et al. 2012;Robbins and Monro 1951) or ADAM (Kingma and Ba 2015).In any case, the solution to * is single valued, and any new prediction made by the trained NN at a new input ŷnew = f (x new ; * ) would also be single valued.No uncertainty information is produced to convey the confidence or quality of these estimates that would be affected by, for example, how many points were used to train the NN and how noisy those training points were.Having uncertainty quantification for the model and its predictions will be particular important to support safety-and missioncritical decision-making, such as in icing detection.

Bayesian neural networks
We adopt a probabilistic approach that seeks to produce not just a single-valued prediction, but a distribution of plausible prediction values.This distinction is demonstrated in Fig. 3.In particular, we take a Bayesian approach and treat as a random variable with an associated probability density function (PDF) that reflects the uncertainty  where p( | x T , y T ) is the posterior PDF on the weight parameters; p( ) is the prior PDF; p(y T | x T , ) is the likeli- hood PDF; and p(y T | x T ) is a normalizing term independent of often known as the model evidence or marginal likelihood.Solving the Bayesian inference problem then constitutes computing or characterizing the posterior p( | x T , y T ) for a given training dataset (x T , y T ) .Once the posterior is available, we can propagate its uncertainty to predictions via Monte Carlo sampling (i) ∼ p( | x T , y T ) and correspond- ingly produce samples of ŷ(i) new = f (x new ; (i) ) that represent the posterior predictive uncertainty-that is, a distribution of plausible prediction values-at any new input x new .
Solving for the posterior is extremely difficult for NNs, since can easily reach thousands, even millions of dimension for large NNs intended to model complex scientific and engineering processes.While Markov chain Monte Carlo algorithms (Gilks et al. 1996;Andrieu et al. 2003;Robert and Casella 2004;Brooks et al. 2011) are typically used for Bayesian computations for their ability to generate samples from the exact posterior, these methods do not scale well to such high dimensions in practice.In contrast, variational inference (VI) (Blei et al. 2017;Zhang et al. 2019) formulates an optimization problem that seeks the best approximation to the posterior from within a class of parameterized distributions in lieu of sampling the posterior.If we use q( ; ) to denote the approximate posterior with representing the parameterization of these approximating distributions (e.g., if the q is from the family of all Gaussian distributions, then is the mean and covariance Σ of a Gaussian) and the VI formulation then seeks * that minimizes the Kullback-Leibler divergence (a measure of farness between two distributions) from the true posterior to the approximate posterior: Since an optimization task is only concerned with finding the single optimum point, it is more scalable to higher dimensions than a sampling task that needs to portray the entire landscape of the posterior distribution.Hence, by converting a sampling task to an optimization one, VI effectively trades off accuracy of posterior representation for scalability, making it more suitable for BNNs.

Projected stein variational gradient descent
We focus on two particle-based VI methods: Stein variational gradient descent (SVGD) (Liu et al. 2016) and projected SVGD (pSVGD) (Chen and Ghattas 2020).With particle representations, these methods are able to capture correlation and non-Gaussian structures of distributions that the commonly used mean-field Gaussian VI cannot.
SVGD is derived by leveraging the relationship between the (functional) gradient of objective in Eq. ( 5) to the Stein discrepancy, the latter which can be approximated using a set of particles (see Liu et al. 2016 for a detailed derivation).A gradient descent procedure can then be formed to iteratively minimize this gradient, as summarized by where i , i = 1, … , m is the location of the ith particle at optimization iteration , is the learning rate that can also be adapted with , and φ * () is the update direction.In turn, the update direction can be expanded as with k(⋅, ⋅) being a positive definite kernel.Furthermore, the gradient of the true log-posterior in the above equation can be evaluated via the sum of gradients of log-likelihood and log-prior, since the gradient of the log model evidence with respect to is zero.The overall effect is an iterative transport of a set of particles to best match the target posterior distribution p( | x T , y T ) .One observed drawback to SVGD, however, is the tendency for particles to collapse toward the mode of the distribution as the dimensionality becomes high relative to the number of SVGD particles, thereby (potentially drastically) under-representing the uncertainty (Chen and Ghattas 2020).
pSVGD (Chen and Ghattas 2020) aims to address this issue by projecting parameters into a low-dimensional subspace and conducting the particle-based inference in this subspace.This allows inference to take place with a much higher particle-to-dimension ratio, which can mitigate the posterior collapsing phenomenon from SVGD.The specific space selected for projection is a likelihoodinformed subspace (LIS).Components with the greatest likelihood-to-prior ratio are identified for inference while components with a negligible likelihood-to-prior ratio are assumed to be unmodified from the prior and need not be inferred directly.In other words, we will focus and work in only the subspace where the data can make a significant update from prior to posterior.We refer readers to Chen (6) Page 7 of 17 197 and Ghattas (2020) for a detailed derivation and present a brief description below.
The LIS is determined by first calculating the average (log)-likelihood Hessian − ∇ 2 L where L = log p(x T , y T | ) .−∇ 2 L can be calculated either directly via back-propagation of the NN to evaluate the second derivative or much more inexpensively via the Gauss-Newton approximation −∇ 2 L ≈ −(∇ L∇ L T ) needing only first- derivative information.Next, we solve the generalized eigendecomposition problem: where i are the r leading eigenvalues that represent the relative influence of the likelihood over the prior, i are the corresponding eigenvectors, and Γ −1 0 is the inverse of the prior covariance.Therefore, eigenvalues below a certain threshold correspond to components with negligible influence from the prior, leaving the prior effectively unchanged.These negligible eigenmodes can then be truncated, and the LIS is constructed with the surviving eigenvectors corresponding to non-negligible eigenvalues.The choice of threshold value for truncation is application specific.For lower dimensional, mostly linear models, a threshold value O(1) is sufficient.For larger, more complex models, a lower threshold value may be necessary.The original authors of pSVGD (Chen and Ghattas 2020) use a value of 0.01.
With the LIS found, the prior parameters 0 can be decomposed into the projected component r 0 and the its complement T 0 (subscript 0 indicates that the decomposition is done under the prior distribution).The complement T 0 is then frozen at the prior projection.The projected component of the parameters will be inferred under the low-dimensional projected variable using the same iterative procedure in SVGD, where r = Ψ r .The original full-parameter vari- able after pSVGD inference can then be reconstructed at the end: Performing Bayesian inference on instead of permits a much lower-dimensional problem.There is also secondary benefit that by truncating to the r components most informed by data, the complement components that would otherwise be more amenable to fitting to data noise are now regularized, further improving predictive performance.

Ice prediction
The PoliMIce (Gori et al. 2015a) 4 are in relatively good agreement with the experimental measurements.The ice predictions at the lower temperature conditions displayed in Fig. 4a and b closely resemble the ice measurements.There are slightly greater discrepancies when moving to the higher temperatures as shown in Fig. 4c and d.The numerical predictions tend to predict the presence of glaze ice earlier than the experiments.Proceeding the icing validation assessment, the subsequent goal of this work was to generate a dataset of ice shapes for the NN model.The experimental work from Shin and Bond was selected with this in mind.In their work they studied the effect of 7 different temperatures on the ice shape.Due to the NN model requiring a significantly larger dataset, this work increments the temperature by 0.5K to provide 50 samples.Additionally, ice shapes at a MVD of 40 μm are simulated to provide 100 samples in total.The icing conditions used for all the samples are summarized in Table 2.

Mesh generation
The focus on modeling small-scale vortical structures using scale-resolving simulations capabilities in SU2 requires three-dimensional grids.An in-house grid generation software called uhMesh (Dussin et al. 2009) is utilized to automatically generate a three-dimensional mesh from a two-dimensional ice shape.This is achieved by extruding the ice shape for 1 chord length in the spanwise direction.
The assumption is that while the ice geometry contains no variation in the spanwise direction, the flow does due the scale-resolving modeling approach.Simulating a fully threedimensional icing dataset of this size would be beyond the current hardware capabilities.Additionally, the ice prediction approach based on Myers' model (Myers 2001) would not capture the spanwise variation in the ice shape.To capture such three-dimensional roughness in the ice shape would require an approach, such as the one suggested by Szilder and Lozowski (2018).However this would be highly challenging to mesh in order to determine the aerodynamics and aeroacoustics of the iced wing.
Figure 5 illustrates the mesh generation for the clean wing using this approach.With scale-resolving simulations in mind, a boundary layer with a height of approximately 0.5c is used to generate a structured region surrounding the wing.The boundary layer at the wall is resolved to ensure y + < 1 everywhere.The mesh is extruded in the spanwise direction for 1 chord length and discretized by 150 elements.The total number of elements for each mesh is approximately 16 M with it varying by as much as 500k depending on the ice shape.

Turbulent flow simulation and noise prediction for different ice shapes
To compare the turbulent flow field and far-field noise features of the clean and iced wings, the clean NACA0012 airfoil as well as a rime-and a glaze-iced airfoils associated with Rime 1 and Glaze 2 in Table 1 are presented.The scale-resolving EDDES-FWH approach outlined in Sect.colored by non-dimensionalized velocity magnitude, as shown in Fig. 6.The effect of leading-edge icing is evident.
While the clean and rime ice airfoil show mostly attached flows over the suction side, detached eddies signifying separated flow can clearly be seen in the more severely iced (glaze ice) airfoil, where the Q-criterion shows more fine-scale, three-dimensional turbulent structures over the entire suction side.Note that even though the ice shapes are quasi-2D, the resolved turbulent structures by EDDES are very much three-dimensional, especially in the glaze ice case where the Kelvin-Helmholtz instability is triggered immediately past the ice horn and develops quickly into 3D turbulence.
The time histories of the lift, drag, and moment coefficients of the rime and glaze ice wings, computed using the EDDES-SA simulations are compared in Fig. 7.In the case of rime ice, the force and moment coefficients are very similar to those of the clean airfoil.Small differences between the clean and rime ice only occur at the peaks of the values which are associated with the maximum and minimum angle of attack.On the contrary, the glaze ice airfoil exhibits significantly deteriorated lift and strongly elevated drag, signifying flow separation over a large portion of the airfoil suction side.The moment coefficient is also severely effected during the glaze ice regime.This is confirmed by examining the wall shear stress of the three cases shown in Fig. 8.Both clean and rime ice airfoils show trailing edge separation due to the dynamic effect of the upstroke.Surface noise footprints visualized by the instantaneous wall pressure fluctuation ( p � = p − p ) of the clean, rime, and glaze ice airfoil sections are compared in Fig. 9.The clean and rime ice airfoils exhibit similar surface noise footprints.On the other hand, the glaze airfoil exhibits a highly irregular noise footprint pattern  with significant spatial frequency content (Fig. 10).One may already speculate that the far-field noise signature of the glaze ice should be associated with elevated highfrequency contents.This is confirmed by the far-field noise spectra shown in Fig. 11-while the rime ice airfoil shows a slight increase of high-frequency noise over the clean airfoil, the glaze airfoil has a significantly elevated noise level over the entire high-frequency range up to the cut-off.
Clearly, the more severe the leading-edge ice accretion is, the more elevated the high-frequency range of the far-field noise tends to be.This can be exploited as a distinguishing feature in a machine learning model to build a mapping between the far-field noise level and ice shape as well as the associated aerodynamic performance of the ice airfoil.

Machine learning predictions of aerodynamic performance under icing conditions
Following Sect.2.3, we build BNNs to predict aerodynamic performance metrics from inputs of acoustic measurements.Specifically, the BNN takes input of a compressed principal component analysis (PCA) representation of the power spectral density (PSD) of the acoustics time series and predict the minimum, mean, and maximum values of the lift, drag, and moment coefficients over the duration of the time series window.To achieve this, we pre-process the acoustic time series in the training dataset by computing the log values of PSD for each data signal covering the frequency range 17-1000 Hz and perform a PCA to find and truncate to only the 6 leading components.All PSDs are then projected into this lower-dimensional representation.All inputs are then standardized using the mean and standard deviation statistics from across all training data simulations.For the output variables, we extract from each simulation's aerodynamic time series the minimum, maximum, and mean over the time window for C L , C D , and C M , a total of 9 quantities of interest (QoIs).Overall, the BNN maps from a 6-dimensional input to a 9-dimensional output.The precise architecture is shown in Table 3, with one hidden layer of 15 neurons leading to a total of 249 model parameters.We purposefully favor a small BNN architecture due to the relatively few training samples available.
To set up the BNN, the prior p(w) from Eq. ( 4) is chosen to be independent Gaussian w k ∼ N(0, 0.2 2 ) , and the likeli- hood p(y T | x T , w) is chosen to reflect an additive independ- ent Gaussian data noise: We use 600 particles to perform SVGD and pSVGD for our 249-parameter BNN.In order to obtain a robust measure of predictive performance of the BNN, we perform fivefold testing as follows.We divide the overall dataset of 93 simulations into fivefolds of 18-19 simulations each.Then, 5 BNN models are built where each model is trained on 4 of these folds and tested on the remainder fold.This allows us to obtain test performance statistics where each of the 93 data points gets a chance to play the role of a testing point.
Figure 12 presents both the SVGD and pSVGD BNN predictions, with each panel showing an aggregated result from the 5 testing scenarios.Using BNNs, we now have uncertainty information along with predictions.These uncertainties are depicted as error bars (± 1 standard deviation) in the plot, which combine both the epistemic uncertainty calculated via SVGD/pSVGD and the output aleatoric uncertainty specified in the likelihood above (i.e., the additive Gaussian data noise k with = 0.2 ).Note that input aleatoric uncertainty (i.e., measurement error of the acoustic signal) is not considered in this work, and requires further investigation since they are often not additive independent Gaussian.The SVGD and pSVGD results match well overall.pSVGD recovers slightly greater uncertainty of around 16% across all predictions compared to SVGD which is amenable to underpredicting uncertainty due to particle collapse.The reduced dimension from pSVGD is about 80% lower than that of SVGD.While pSVGD uses additional approximation in its eigenvalue truncation, we believe this truncation error is outweighed by the benefit from a greater particleto-dimension ratio resulting from the dimension reduction, leading to an overall more accurate posterior representation for pSVGD.
Posterior predictive expectation of the absolute and normalized root-mean-square error (RMSE) (i.e., posterior-averaged RMSE) for each QoI by the pSVGD model are shown in Table 4, showing good agreement with the true target values.However, we draw caution to interpreting these results, since BNNs are not meant to minimize error directly (error is not its loss function) but rather to accurately capture the uncertainty in the NN as justified by the available training data at hand.As such, while it is still valuable to understand the error performance from these BNNs, one must be careful to not judge these models based on error alone.
In general, we also observe rime ice predictions to achieve lower error and lower uncertainty than glaze ice predictions, which matches our intuition where generally more complex aerodynamic flows are caused by the more severe glaze icing conditions.These trends are shown in

Conclusion
This work introduces an original methodology for developing a Bayesian neural network acoustic ice detection system.The acoustic signatures generated by different iced wings are harnessed to predict their aerodynamic performance.An inter-disciplinary methodology is adopted to develop the Bayesian neural network acoustic ice detection system.Fig. 12 Predicted versus true values for all QoIs, using SVGD (top row) and pSVGD (bottom row).The two sets of results overall match well with each other.More specifically, pSVGD recovered approximately 16% more of the Bayesian predictive uncertainty relative to its SVGD counterpart that is amenable to underpredicting uncertainty due to particle collapse.pSVGD also achieved 80% dimension reduction of the parameter space compared to SVGD The PoliMIce toolkit is utilized to predict a test matrix of iced airfoils.The aerodynamics and aeroacoustics associated with the iced airfoils are, respectively, solved in SU2 using high-fidelity EDDES scale-resolving simulations and the FWH-F1A approach.A Bayesian neural network is then trained using projected Stein variational gradient descent to provide (1) mapping from observed far-field broadband noise levels to predictions of aerodynamic performance indicators and (2) uncertainty information reflecting the neural network's confidence in its predictions.While the icing simulation, aeroacoustic prediction as well as the construction of Bayesian neural network model all require potentially extensive computational resources in the offline phase, the prediction of aerodynamic performance indicators based on measured far-field noise level can be obtained rapidly in the time-critical online phase in-flight.
To the best of our knowledge, this work is the first to employ high-fidelity aeroacoustic simulations to develop correlations between aerodynamic performance and farfield noise level of iced airfoil sections via machine learning techniques.That being said, it should also be noted that the current work as-is only serves as a proof-of-concept study.It is clear that the pitching wing considered in this study is a highly simplified representation of the complex rotorcraft icing problem.Various noise sources such as those from tail rotors and engines, as well as rotor-fuselage interaction must be taken into account and differentiated in developing an ice detection system using the multi-disciplinary framework proposed in the current work.Acknowledgements The authors gratefully acknowledge the highperformance computing resources provided by the UK National Supercomputing Service via ARCHER2 (http:// www.arche r2.ac.uk), on Project e734 (High-Fidelity Aerodynamic and Aeroacoustic Simulations for Urban Air Mobility (UAM) Concepts).JH and XH are supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under Award Number DE-SC0021397.This paper was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Declarations
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.
Responsible Editor: Chao Hu Jeremiah Hauth, Alberto Guardone, Xun Huan, and Beckett Y. Zhou have contributed equally to this work.Topical Collection: Advanced OptimizationEnabling Digital Twin Technology Guest Editors: C Hu, VA González, T Kim, O San, and Z Hu 197 Page 2 of 17

Fig. 1
Fig.1Workflow routes illustrating the physics causation (black solid), a traditional inverse problem route that needs to first identify the ice profile from acoustic signal and then solve for performance indicators (dashed red), and proposed ML model "shortcut" mapping (dotted-dash blue).(Color figure onlilne)

Fig. 2
Fig. 2 Multi-step ice accretion simulation a) Neural network (b) Bayesian neural network

Fig. 3
Fig. 3 Illustration of a simple NN and BNN, each with an input layer, two hidden layers, and an output layer.Trainable parameters are depicted in red as deterministic values in the NN and as probability distributions in the BNN

Fig. 4
Fig.4Comparison of ice predictions with measurements at different temperatures.The ice shape is displayed in non-dimensional form where the coordinates, x and y, are normalized by the wing chord length, c 2.2 is used to model the fine-scale turbulent structures induced by the ice accretion.With ice detection on rotorcraft in mind, the 1c section of the wing is subsequently set to pitch around its quarter-chord to represent the cyclic pitching motion of a rotor blade in forward flight.While this is a primitive way to model a rotor blade in forward flight, it was a requirement due to rotorcraft icing simulations requiring fully three-dimensional simulations which are computationally very costly.The sinusoidal pitching motion of the wing can be described by = 4 • ± 4 • and pitching with a frequency of 6.0 Hz.The mean angle of attack is consistent with the icing experiments presented in Sect.3.1.A time step of Δt = 0.05c∕U ∞ = 0.00025 s is used, which results in 666 time steps within each pitching cycle.The flow statistics and acoustics are computed after the 5th pitching cycle, over 10 pitching cycles.The instantaneous flow fields at AoA = 4.0 • on the upstroke are visualized by the iso-surfaces of the Q-criterion (a) Close-up of the mesh surrounding the wing.(b) Isometric view of the surface mesh.

Fig. 6
Fig. 6 Q-criterion iso-surface for the clean and iced wings colored by non-dimensionalized velocity magnitude, U∕U ∞

Fig. 7
Fig. 7 Time histories of lift, drag, and moment coefficients of the clean, rime, and glaze ice airfoils computing using the EDDES-SA simulations

Fig. 8
Fig. 8 Skin friction coefficient, C F x , for the clean and iced wings

Fig. 9 Fig. 10
Fig. 9 Instantaneous pressure fluctuations, p ′ , for the clean and iced wings

Fig. 13
Fig. 13 MSE of the posterior predictive distributions (left) and standard deviation of the posterior predictive distributions (right) averaged over all QoIs.The predicted QoIs are normalized to allow direct comparison and averaging

Table 2
Test matrix to produce sample of rime and glaze ice shapes

Table 4
The posterior predictive expectations of the normalized and absolute RMSE for the pSVGD model (i.e., posterior-averaged RMSE) Normalized errors are calculated based on transforming the true test target values and predictions by the mean and standard deviation of the true training target values, which in turn can be interpreted roughly as percent error