Machine learning-enabled prediction of wind turbine energy yield losses due to general blade leading edge erosion

Blade leading edge erosion is acknowledged to significantly reduce the energy yield of wind turbines. The problem is particularly severe for offshore installations, due to faster erosion progression boosted by harsh environmental conditions. This study presents and demonstrates an experimentally validated simulation-based technology for rapidly and accurately estimating wind turbine energy yield losses due to general leading edge erosion. The technology combines the predictive accuracy of twoand three-dimensional Navier–Stokes computational fluid dynamics with the runtime reductions enabled by artificial neural networks and wind turbine engineering codes using the blade element momentum theory. The main demonstration is based on the assessment of the annual energy yield of the National Renewable Energy Laboratory 5 MW reference turbine affected by leading edge erosion damage of increasing severity, considering damages based on available laser scans and previous leading edge erosion analysis. Results also include sensitivity studies of the energy loss to the wind characteristics of the installation site. It is found that the annual energy loss varies between about 0.3 and 4%, depending on the damage severity and the site wind characteristics. The study also illustrates the necessity of resolving the geometry of eroded leading edges rather than accounting for the effects of erosion with surrogate models, since, after an initial increase of distributed surface roughness, erosion yields leading edge geometry alterations causing aerodynamic losses exceeding those due to the loss of boundary layer laminarity consequent to roughness-induced transition. The presented technology enables estimating in a few minutes the amount of energy lost to erosion for many-turbine wind farms, and offers a key tool for predictive maintenance.


Introduction
Blade leading edge erosion (LEE) of wind turbines is an important and unresolved problem in the wind energy industry. Blade LEE spoils the aerodynamic performance of wind turbine rotors, causing reduced turbine power and, ultimately, losses of annual energy production (AEP). The maximum tip speed of modern offshore wind turbine blades reaches values of 90 to 100 m/s, and higher tip speeds may help further reduce the levelized cost of energy of wind-based electricity, due primarily to lower cost of gearboxes and drivetrains. This cost reduction is enabled by the reduced torques on the low-speed shaft required to transmit a given power amount [1]. At the high velocities of the outboard part of the blade, the leading edge (LE) is progressively eroded due its collision with rain, hail and sand. The damage alters the geometry of the blade airfoils and spoils their aerodynamic performance, resulting in increased drag and reduced lift forces of the affected blade sections.
Several LE protection solutions have been proposed by the industry and academia [2][3][4][5], and considerable research is presently ongoing in this area worldwide. Although in some cases there is some evidence that these solutions may mitigate LEE, it is still unclear if any of these can reduce below acceptable thresholds the detrimental impact of LEE on wind turbine and farm energy yield. The resulting uncertainty makes it impossible to adopt a predictive maintenance approach for LEE, which remains the largest cause of unscheduled blade maintenance. One key problem hindering the adoption of a predictive maintenance approach stems from the difficulty of estimating reliably the turbine power and AEP losses associated with specific LEE pattern and extent. The availability of this quantitative relationship would give a strong support to scheduling blade surface maintenance and repair. This would enable one to optimize blade maintenance by considering the trade-off of maintenance costs and revenue increase due to not incurring LEEinduced AEP losses. For lift-based machines, such as horizontal axis and Darrieus wind turbines, estimating AEP losses requires knowledge of the power curve of the turbine affected by LEE. Since the power curve depends directly on the lift and drag forces of the blade airfoils, reliably predicting AEP losses requires reliably predicting the lift reduction and the drag increase experienced by the affected blade sections.
A different aerodynamic loss mechanism and magnitude is expected depending on the surface damage pattern and extent, with both factors depending on the stage at which the erosion is. Hence, a detailed description of the damage is required for a reliable estimate of its aerodynamic effect. The generation and progression of LEE are a sitedependent process observed by wind farm operators and reported in several studies [2,[6][7][8]. At the first stage, erosion alters the surface geometry on a microscale, leading to increased surface roughness. When the material reaches the damage incubation limit [9], small pits appear on the surface. At a more advanced stage, the pit size increases, and groups of pits coalesce in larger gouges. At the most critical stage, erosion removes larger amounts of material, along the adhesion boundaries between different material layers and LE delamination occurs. The four subplots of Fig. 1, reported in [10], provides an example of the aforementioned LEE stages using photographs of real blade LEs. Photographs (A) to (D) illustrate the progression from increased surface roughness to advanced LE delamination.
Gaudern et al. [11] proposed a classification of LEE stages based on a campaign of field observations, and studied LEE using wind tunnel testing. Each stage was modeled by applying rough films to the LE of an airfoil model, and lift and drag over a representative range of angle of attack (AoA) were measured to assess the impact of each pattern on the airfoil performance. The effect of LE roughness on blade airfoil performance was investigated using both wind tunnel testing and numerical simulations also in [12,13]. Sareen et al. [14] carried out a comprehensive wind tunnel test campaign considering the DU 96-W-180 wind turbine airfoil subjected to different patterns and severity of geometry alterations. The studied erosion patterns were generated on the basis of field observations and photos. Each LEE stage analyzed in the experiment was characterized by a different distribution of pits, gouges and delamination patterns. Wang et al. [15] studied the effect of pitting erosion using two-dimensional (2D) computational fluid dynamics (CFD) analyses, describing the pits as semicircular cavities distributed around the airfoil LE. Castorrini et al. [16] presented and demonstrated a three-dimensional (3D) CFD method for analyzing the impact of LE pits and gouges on airfoil performance. A discussion on the aerodynamics of airfoils featuring LE erosion cavities, supported also by a validation study based on the measured data of [14], was reported in [17]. The study also highlighted the importance of both modeling smaller-scale surface roughness and resolving geometrically sparse and larger LE erosion cavities when analyzing initial and intermediate LEE stages. Schramm et al. [18] used 2D CFD to assess the impact of a particular LE delamination pattern on the AEP of the National Renewable Energy Laboratory (NREL) 5 MW reference turbine [19]. All these investigations provide a strong foundation for understanding key aspects of the aerodynamics of blades affected by LEE, but do not focus on industrial application, since they refer to specific and/or idealized erosion damages, whereas the LEE-induced power reduction and AEP loss of wind turbines and farms depend on significantly more complex and heterogeneous LEE patterns. The objectives of this study are to a) define a general methodology to estimate turbine power reduction and AEP losses due to general LEE geometry, characterized by pits, gouges and delamination of fairly arbitrary pattern and size, b) demonstrate this methodology to estimate the growth of the AEP loss corresponding to increasing severity of the erosion stage, and c) highlight the necessity of resolving the damage geometry of most erosion stages to avoid underestimating wind turbine and farm AEP losses. The baseline concept of the AEP loss prediction system (ALPS) was presented and demonstrated for a single LEE stage, namely that of severe delamination, in [20]. In that study, ALPS used machine learning (ML) trained with a pre-computed large database of delaminated airfoil force data for the rapid calculation of airfoils with general delamination. The study used a 2D damage geometry parametrization and 2D CFD, as the geometric representation of this damage based on 2D chordwise blade strips is deemed to be a good trade-off of computational costs and modeling realism. In the present study, the ALPS technology is generalized to enable the analysis of additional LEE stages, namely those of initial and intermediate LEE, characterized by pits and gouges. These cavities are inherently 3D geometric features whose CFD analysis requires more sophisticated 3D geometry parametrization and mesh generation, and computationally more costly 3D CFD simulations. To demonstrate the application of the generalized ALPS technology, the method is applied to determine the power reduction and AEP loss of a utility-scale wind turbine caused by a set of general LEE damages of increasing severity. This enables estimating the growth of AEP losses during operation before blade surface maintenance takes place. The LEE analyses reported herein focus particularly on the initial/ intermediate LEE stages, characterized by increased surface roughness first, and pits and gouges successively (top images of Fig. 1), and more advanced LEE, characterized by severe LE delamination (bottom right image of Fig. 1).
Section 2 summarizes the ALPS concept, while Section 3 defines the CFD and ML methods for the prediction of lift and drag coefficients of blade sections featuring LE pits, gouges and delamination. The automated geometry parametrization of the 3D airfoil surface with LE erosion cavities and the CFD grid generation are reported in Section 4, which also provides the geometric definition of all pitted and gouged airfoils included in part of the ALPS force databases of this study. Section 5 has three parts: the first presents the validation of the adopted CFD method against recently published measured force data of nominal and delaminated airfoil geometries; the second part discusses with theory and examples the sensitivity of the laminar-to-turbulent transition of the blade airfoil boundary layers (BLs) to LE erosion features, and motivates the choice of the fully turbulent flow assumption made to generate the ALPS force databases; the third part is a validation study of the ALPS ML predictive capabilities. The key results of this report are in Section 6: the first part provides the geometric definition of specific LEE damages, ordered by severity; in the second part ALPS is used to estimate the progressive degradation of the power curve of the NREL 5 MW turbine subjected to these LEE damages; the third part is a parametric analysis of the dependence of the AEP loss on both LEE damage type and extent, and wind characteristics of the installation site; and the fourth part is an sensitivity analysis quantifying the uncertainty on AEP losses due to uncertainty in the estimate of the airfoil force coefficients. A summary of the study and future work are provided in the closing Section 7.

AEP loss prediction system
This section, summarizes the ALPS analysis framework [20,21], with particular emphasis on its main components and dependencies. The ALPS modular structure is provided in graphical form by the blockdiagram of Fig. 2, which is used to support the system description below.
The main purpose of ALPS is to determine the AEP of a wind turbine whose blades are affected by LE erosion damage of general shape and severity (block 1 of Fig. 2). The estimate of the damaged turbine AEP requires determining its power curve first. The power curve can be determined using a wind turbine engineering code using the blade element momentum theory (BEMT) for blade aerodynamics. The BEMTbased calculation of the power curve of the damaged turbine requires knowledge of the lift coefficient c l and drag coefficient c d of the blade sections over the whole range of relative AoA encountered along the blade length. The calculation of c l and c d starts from large databases of aerodynamic force coefficients of damaged airfoils (block 2), precomputed with 2D and 3D CFD in the present ALPS implementation. One database is generated for each damage class, e.g. one for pits, one for gouges and one for delamination. Each damaged airfoil of a database features c l and c d curves defined over its typical AoA operating range. In general, the databases will not contain the eroded airfoil geometries of the case at hand, but rather airfoils with erosion patterns and magnitudes of the same class. Therefore, ML is used to exploit the database knowledge to rapidly get the c l and c d data of the damaged blade airfoils at hand for the required range of AoA (block 3), with no need of further CFD analyses. Artificial neural networks (ANNs) trained and tested using the databases are presently used in ALPS, as discussed in further detail in Section 3. The force coefficients of the eroded blade airfoils at hand are then used by a wind turbine aerodynamic analysis code based on BEMT or a multi-disciplinary engineering code, such as the NREL OpenFAST code [22], using BEMT for aerodynamics (block 4). Using a given power and load control algorithm (block 5), the wind turbine code determines the wind turbine power and loads over the entire range of operating wind speeds (block 6). Finally, the AEP estimator (block 7) integrates the turbine power curve against the annual wind frequency distribution at the rotor hub associated with installation site (block 8) to determine the sought AEP of the damaged turbine (block 9).
The current implementation of ALPS also requires as input the nominal geometry of the turbine blade and some information about the turbine control, such as minimum and maximum angular speeds of the rotor. The construction of the databases of damaged airfoils requires selecting suitable range limits of all geometric parameters defining the LE erosion damages under investigation. For the case of LE delamination examined in [20], the ranges of all delamination parameters in the database were selected on the basis of damage geometry information available in [14] and other published data. The ranges of all geometric parameters defining the LE damage by erosion cavities in [16] were determined using the findings of a survey of eroded LE photographs available in the published literature and data reported by the industry. The specific LE damages by erosion cavities and delamination examined in this study (Section 6) using ALPS are generated using a more general and comprehensive approach, reflecting also the alterations of the LE geometry as LEE progresses during operation. This approach is based both on surveys of available LEE photographs and laser scans, and numerical estimations of LE material loss due to erosion [23][24][25].
In ALPS, the turbine control of the nominal and damaged turbines can be set to be the same. However, the alterations of the aerodynamic characteristics of eroded blades may also result in alterations of the optimal tip-speed-ratio (TSR), namely the value of this parameter at which the maximum value of the turbine power is achieved. Because of this phenomenon, further discussed in [16], an adaptive power control aiming at mitigating power losses due to LEE was proposed in [20]. The method tracks the variations of the optimal TSR due to LEE and enables the turbine with LEE to operate at the maximum possible power before the rated wind speed is achieved. Above the rated wind speed, LEE does not spoil the turbine power, as the rated power can be achieved by altering the AoA using the blade pitch control, as discussed in [20].

Numerical method
The main computational technologies used in this study are computational fluid dynamics and machine learning. The specific algorithms and codes are summarized in the two following subsections.

Navier-Stokes computational fluid dynamics
The lift and force data past airfoils with LEE are determined by simulating their turbulent viscous flow field. The Navier-Stokes CFD code used for all 2D and 3D aerodynamics simulations reported herein is version 2019 -release 3 of ANSYS FLUENT [26]. The air flow is modeled as incompressible and single-phase. The pressure-based steady Reynolds-averaged Navier-Stokes (RANS) equations are solved in all cases, and the turbulence closure is accomplished by means of Menter's two-equation k − ω shear stress transport (SST) model. [27]. The space discretization of all conservation laws uses the second order upwind scheme. The COUPLED solver, whereby the continuity and momentum equations are solved in a strongly coupled fashion and all other transport equations are solved in a loosely coupled fashion, is adopted for the integration of the governing equations.
Modeling of the laminar-to-turbulent transition of the blade airfoil BLs relies on the use of two additional transport equations coupled to the SST turbulence model. The resulting model is the Langtry-Menter fourequation transitional SST model [28][29][30], referred to as the γ − Re θ SST model hereafter. The method couples the SST turbulence model with socalled local correlation-based transition modeling. This approach does not simulate the physics of laminar-to-turbulent transition; rather it relies on empirical correlations that relate local BL parameters, such as the momentum thickness θ and other local properties such as pressure gradient and turbulence intermittency γ, to position the transition onset. The wall and far field boundary conditions for all four equations of the γ − Re θ SST model are described in [30].

Machine learning
Three classes of LE damage by erosion are considered in this study, namely small cavities or pits, large cavities or gouges, and delamination. For each of the typically 6 to 8 nominal airfoils of a wind turbine and for each damage class, a large number of CFD simulations of eroded airfoils is carried out to generate databases of aerodynamic force coefficients, whereby output values of the lift coefficient c l and drag coefficient c d are associated with input values of the geometric damage parameters and the angle of attack. The database is then used to train ML algorithmsartificial neural networks in the present study -and generate the regression to correlate input and output variables. The availability of this regression eliminates the need for new CFD simulations to estimate the force coefficients of erosion damages not included in the database.
For all considered damage classes, generic multi-layer perception (MLP) feed-forward ANNs with one hidden layer are used to learn from the c l and c d data patterns in the database, and infer the aerodynamic forces of eroded airfoils not contained in the database. Two separate ANNs for predicting c l and c d are generated for each nominal airfoil and damage class. This approach follows that presented in [20].
One-layer MLPs consist of universal function approximators f g (x) : R n d →R no [31], where n d is the size of the array x of input variables, and n o is the size of the array f g of functional outputs. In the application of this study, the elements of x are the angle of attack of the airfoil oncoming flow and geometric parameters defining the eroded airfoil geometry, as explained in Section 4, and n d is thus the number of parameters defining a particular LE damage geometry. In all cases, n o = 1, since the output function is scalar, corresponding to either c l or c d .
The general matrix-vector definition of f g is: where W (1) is a weight matrix of size (N × n d ) and N is the number of neurons in the hidden layer, b (1) is a bias term represented by a column vector of length N, W (2) is a weight matrix of size (n o × N), b (2) is a bias term represented by a column vector of length n o , and A1 and A2 are the activation functions of the hidden layer and the output layer, respectively. Typically, A1 is the hyperbolic tangent sigmoid function, defined as tanh(a) = (e a − e − a )/(e a + e − a ), and A2 is linear. With this choice, function f g (x) becomes: The hyperbolic tangent sigmoid function is a well established activation function in MLP feed-forward ANNs for data fitting, and, in light of the good results obtained with this choice in the presented analyses, no alternatives have been investigated.
Given a set of N s training samples {(x 1 , y 1 ), …, (x Ns , y Ns )}, where y i is the observed response to the input x i , a learning algorithm seeks the values of the weight matrices and bias vectors that minimize the difference between part or all of the observed N s sample responses y i and the N s responses given by f g (x i ). Gradient-based optimization is used to determine the weight matrices and the bias vectors by minimizing the aforementioned difference. Since the responses are assumed to be smooth functions of the inputs x i and the internal weights and biases, the gradients of the difference between sample responses and ML approximations with respect to weights and biases can be computed using the so-called back-propagation method, which relies on applying the chain rule for derivation. The back-propagation approach is applied repeatedly to propagate gradients through all layers, starting from the output at the top, where the network gives its response, and working all the way down to the hidden layer at the bottom, where the input is provided. Once the gradients with respect to weights and biases for each layer are computed, the objective function expressing the level of fitting of the training data can be optimized.
If the ANN has a sufficiently high number of degrees of freedom, a parameter determined by the selected number of neurons for the hidden layer, the fitting error can be reduced to machine zero; this results in the system overfitting the training data and being possibly unable to generalize adequately its predictions to regions of the input space where there are insufficient or no training data. To mitigate this risk, two different approaches have been adopted in this study. One approach is to subdivide the available data set into a training set and a validation set. At each step of the gradient-based optimization via back-propagation (training), the fitting error based on the validation set is also computed. The fitting error is taken to be the mean square difference between the observed responses (CFD lift and drag coefficients) and their ML estimates. Training is stopped when the fitting error based on the validation set starts increasing. This typically happens before the fitting error based on the training set achieves machine zero. The alternative approach introduces one or more regularization terms that modify the fitting function so that the obtained ANN has good generalization properties. One of the most widely used regularization approaches is Bayesian regularization, in which a modified linear combination of fitting errors, weights and biases is minimized [32].
Using either approach to prevent overfitting, the learning process via back-propagation consists of a gradient-based optimization problem, which is solved with the Levenberg-Marquardt (LM) method [33]. Because of the local nature of gradient-based optimization, however, the learning outcome can be affected by the process initialization. To reduce this risk, 100 training runs with different random initialization are carried out both with the training and validation set approach, and the Bayesian regularization approach. The selected ANNs are those which provide a good trade-off of high fitting degree of the training set and high generalization properties on the test set. In the parametric analysis leading to the selection of the best suited ANNs, the N parameter is also varied to determine the value that is sufficiently large to ensure good ANN predictive capabilities, but not unnecessarily large, so as to limit the computational cost of the ANN construction. The optimal value of N for both c l and c d ANNs for LE erosion damage by pits and gouges of the present study varies between 5 and 6, whereas that of the ANNs for the case of LE delamination presented in [20] varied between 20 and 30.
For the ANNs predicting force coefficients for damage by LE gouges, the training set is 85% of the randomized database when using Bayesian regularization, and the remaining 15% forms the test set, used to test the generalization strength of the ANNs. The size of the training set is instead 68% of the randomized database and that of the validation set is 17% of the same database when using the approach without regularization, and the remaining 15% forms the test set. For the ANNs predicting force coefficients for damage by LE pits, the size of the training set is 97% of the randomized database when using Bayesian regularization, and the remaining 3% forms the test set. The size of the training set is 78% of the randomized database and that of the validation set is 19% of the same database when using the approach without regularization, and the remaining 3% forms the test set.

Geometry, mesh and database definition
This section presents the numerical set-ups for the aerodynamic performance analysis of the NACA 64-618 airfoil with LE erosion cavities. This airfoil is used from 70% rotor radius to the blade tip of the NREL 5 MW reference turbine, which is used for the AEP analyses of Section 6. The first subsection describes the geometry parametrization and CFD grid generation of the airfoils featuring LE erosion cavities included in the ALPS databases created for this study; the second subsection provides the values of the geometric parameters defining all considered eroded airfoils, and all key parameters used for the generation of the force coefficient ANNs.

Airfoil geometry, grid generation and flow regime
Each eroded section or strip is represented by a 3D airfoil with unitary chord c, small spanwise width, and semispherical cavities in the LE surface with centers on the plane of symmetry of the 3D airfoil. The cavities of a given section have the same radius, and the curvilinear distance of adjacent cavity centers is constant. The choice of the spherical shape for the erosion cavities is made because this geometry was used in the validation of the ALPS CFD set-up against measured force coefficients of an airfoil with this type of LE damage reported in [17]. The parametric assessment of the dependence of aerodynamic performance losses on the cavity shape, also reported in [17], shows a relatively low sensitivity of performance losses to the cavity shape. The geometry parametrization of an eroded LE strip is depicted in the two schematics of Fig. 3, with the left schematic reporting a blade crosssectional view, and the right schematic reporting a planform view. The parameter Φ denotes the diameter of the cavities, whereas the parameters s l and s u denote the curvilinear length of the erosion damage on the lower and upper side of the airfoil, respectively. The latter two parameters indicate the maximum curvilinear length covered by cavities on either side of the LE. The parameter k 1 is the cavity pitch, defined as the curvilinear distance between two adjacent cavity centers normalized by Φ, and is constant for a given strip. The parameter k 2 is the spanwise width of the strip, also normalized by Φ.
The 3D solid geometry of the physical domain for the CFD analysis is generated in three steps. The first consists of defining the 2D spline of the nominal airfoil profile and the complete far field boundary, located at about 30 chords in all directions, in the xy plane. The second step consists of linearly extruding the aforementioned 2D domain in the z direction over a distance k 2 Φ. At the third step, one performs the boolean union of the resulting solid with a set of spheres of radius Φ/2 and centers on the strip symmetry profile, according to the values of s l , s u , and k 1 at hand. A 2D view of the physical domain and an enlarged 3D view of the eroded LE are reported in Figs. 4a and 4b, respectively. The generation of all 3D airfoil models and physical domains is controlled by a single MATLAB script, that reads in all LE erosion data and the nominal airfoil geometry, and, for each damaged airfoil, writes in a suitable format all LE erosion damage parameters, and calls the ANSYS Space-CLAIM CAD preprocessor to generate the boundaries of the entire physical domain. SpaceCLAIM uses a Python script to read in the MATLAB damaged geometry data.
All meshes are generated automatically by using a single MATLAB script to run the generation of all meshes, and a JavaScript file to control mesh settings of the ANSYS MESHING tool. An unstructured grid region discretizes the volume in the neighborhood of the LE cavities, bounded by the cyan surfaces in Fig. 4b. This grid region is a hybrid mesh with a wall inflation layer made up of 25 layers of prismatic elements. The first layer height is 3.9 × 10 − 6 c and the adopted growth rate of the prismatic inflation layer is 1.2. Above the inflation layer, there is a buffer region of the hybrid mesh made up of tetrahedra and pyramids. The latter element type is used at the junction of the buffer region and the remainder of the grid, which extends until the far field boundary and is made up of hexahedral elements. Fig. 4c presents a view of the hexahedral mesh in the symmetry plane of the airfoil strip, whereas an enlarged view of the LE surface mesh in the hybrid mesh region is provided in Fig. 4d.
The geometry parametrization of the LE delamination damage is shown in Fig. 5, and the automated procedure for the damage airfoil geometry construction and grid generation is described in [20].
The number of grid cells of the meshes for airfoils with erosion cavities has been determined on the basis of mesh refinement analyses (not reported for brevity). Mesh-independent c l and c d estimates for values of the AoA α from − 6 • to 12 • have been obtained using meshes with 500 K to 1.5 M cells for about 80% of the considered eroded airfoils. The overall cell number varies significantly with the considered LE damage because the mesh resolution around each cavity and in the spanwise direction is comparable for all geometries, but the number of holes and the model span of each eroded airfoil vary significantly. Meshindependent solutions of the nominal airfoil model with a span of 12 mm have been obtained with a grid of about 500 K cells. For all NACA 64-618 analyses of this article, the freestream turbulence intensity and the eddy viscosity ratio are 5% and 1, respectively, whereas the Reynolds number (Re) is 7 M. The adopted wall distance of 3.9 × 10 − 6 c of  the first grid nodes off the airfoil surface from the airfoil surface itself ensures a maximum nondimensionalized wall distance y + well below 1 around the airfoil for all considered AoA values. Unless otherwise stated, all nominal NACA 64-618 analyses use the γ − Re θ SST model to account for the effects of laminar-to-turbulent transition of the airfoil BLs, whereas all analyses of the eroded NACA 64-618 use fully turbulent BL model, for reasons explained in Section 5.2.

Database definition
A fully automated approach is used to generate the databases of NACA 64-618 airfoils with LE erosion cavities. One database for the pits and one for the gouges are constructed. Each airfoil geometry of a database is defined by the set of 5 geometric parameters Φ/c, s l /c, s u /c, k 1 and k 2 . The value of these parameters defining the airfoils affected by pits and gouges are provided, respectively, in the second and third column of Table 1, and its last row shows that 960 and 192 damaged airfoils are considered for generating the pitted and gouged airfoil databases, respectively. For each damaged airfoil, the CFD analysis used to determine the c l and c d coefficients has not modeled BL transition, for reasons discussed in Section 5.2, and has been run for 10 AoA values, namely from α = − 6 • to α = 12 • with a step of 2 • . Due to numerical instabilities encountered with the flow simulations, however, only 792 pitted airfoils and 179 gouged airfoils with geometry defined by the parameters in Table 1 have been used for the two databases. On the other hand, a test set of 31 gouged airfoils and one with 23 pitted airfoils, all with geometry not included in those of Table 1 have also been analyzed. Therefore, since the input variables on each line of a database include the 5 geometric parameters and also the AoA α, the database of the pitted airfoils has 8150 entries, and that of the gouged airfoils has 2100 entries.
In light of the above, the length n d of the array x of input variables in Eq. (1) is 6. The length n o of the array f g of outputs is 1, since separate ANNs are generated for c l and c d . The number of elements of the training set of the c l and c d ANNs of the pitted airfoils is 7920, and that of the gouged airfoils is 1790. The number of neurons N of the ANNs for predicting the force coefficients of the airfoil with LE cavities varies between 5 and 6.
Using one 16-core node of the HEC computer cluster of Lancaster University [34], one 3D CFD simulation using the grid settings reported above requires about 36 min of wall-clock time, and the construction of the lift and drag curves for each airfoil requires about 6 h of wall-clock time. Thus, the generation of the 1025-airfoil force database requires about 6150 h of wall-clock time using one cluster node. However, the database generation is accomplished by using concurrently 20 cluster nodes, and the required wall-clock time is therefore reduced to about 13 days.
The c l and c d ANNs of the NACA 64-618 airfoil with LE delamination used in the present study are those presented in [20]. LE delamination is parametrized by the three geometric parameters s l /c, s u /c, and d/c, as depicted in Fig. 5, and therefore n d = 4. The number of neurons N of the ANNs for predicting the force coefficients of the airfoil with LE delamination varied between 28 and 30. The bounds of the geometric variables defining the training and validation sets of the delamination database are reported in [20]. The CFD simulations used to generate this database were also fully turbulent, and this choice is also discussed in Section 5.2.

Validation and sensitivity analyses
The first part of this section presents and discusses validation analyses of the CFD methods used to generate the ALPS databases of aerodynamic force coefficients. The second subsection discusses the dependence of the BL state of eroded blade sections (transitional or fully turbulent) on the LE erosion depth and the blade section Reynolds number, with the aim of motivating the choice of the fully turbulent BL model for generating the ALPS databases. The third and final part presents a validation study of the ALPS predictive capabilities for the case of LE erosion cavities, accomplished by comparing the ANN predictions of the aerodynamic force coefficients of a set of eroded airfoils with the values obtained with complete CFD analyses of the same airfoils.

Validation of computational fluid dynamics predictions
The validation of the ALPS CFD set-up for the analysis of the LE damage by erosion cavities using the same mesh generation set-up reported in Section 4 is reported in [17]. That study considers the measured aerodynamic performance of the nominal DU 96-W-180 airfoil and a variant with LE erosion cavities (damage Type B -Stage 2) reported in [14]. The considered flow fields have Re = 1.5 M, and the reported analyses show a good agreement of measured and 3D CFDbased estimates of c l and c d data of these two airfoils (Fig. 6 of [17]). The analyses of [17] also consider the combined effect of distributed surface roughness and greater airfoil geometry damage on the airfoil performance impairment, and highlight the importance of resolving the geometry of the LE erosion damage, rather than accounting for its effect by merely using fully turbulent simulations of the nominal airfoil, also for the case of damage by erosion cavities.
The validation of the ALPS CFD set-up for the analysis of the LE damage by delamination considers the measured aerodynamic performance of the nominal NACA 63 3 -418 airfoil and a variant with LE delamination reported in [35]. The wind tunnel tests have Re = 5 M and turbulence intensity of 0.1%, and the airfoil chord is 1 m. The delamination patch of the damaged airfoil extends from the LE to 3%c on both suction and pressure sides, has height of 0.3%c, and also features distributed roughness of 415 μm. Following the outcome of mesh independence analyses, the grids used for the simulation of the nominal and eroded airfoils discussed below have 156,604 and 172,750 quadrilateral cells, respectively, and have been generated as reported in [20]. The adopted wall distance of 5.  airfoil, the minimum measured and computed LEE-induced c d increases are 0.0072 and 0.0079 respectively, whereas the maximum measured and computed c d increases are 0.116 and 0.0139 respectively. In the same interval, CFD tends to overestimate the lift reduction, since the minimum measured and computed c l reductions are 0.0522 and 0.1144 respectively, whereas the maximum measured and computed c l reductions are 0.1346 and 0.2233 respectively. The impact of these differences on turbine energy yield estimates is assessed in the sensitivity study of Section 6.4. The fully turbulent flow analysis of the delaminated NACA 63 3 -418 airfoil, not reported for brevity, produces c l and c d estimates which differ negligibly from their transitional counterparts. This is because the flow perturbation due to the considered damage produces a flow perturbation sufficient to trigger BL transition at the LE. This phenomenon depends both on the boundary layer thickness and the damage severity. The dependence of LE damage-induced transition on damage depth and freestream Reynolds number is discussed in Section 5.2, which explains why fully turbulent, rather than transitional, flow simulations have been used to generate the ALPS force databases of this study.
The two subplots of Fig. 6 also report c l and c d data of the nominal airfoil computed with a fully turbulent analysis (curves labeled 'nom. ft.'). Comparing these curves with those obtained by simulating the flow of the delaminated airfoil highlights that the performance of the delaminated airfoil is poorer than that of the nominal airfoil with BLs tripped at the LE over the entire AoA range considered. Consequently, analyses assuming fully turbulent BLs of the nominal airfoil to estimate the performance loss due to delamination would underestimate aerodynamic losses. This is because the effect of the forward facing step associated with LE delamination is not only to trip BLs at the LE, but also adding additional profile losses due to thicker BLs [20]. All this illustrates the necessity of resolving the erosion damage geometry, rather than accounting for its effect by merely using fully turbulent simulations of the nominal airfoil, to more reliably estimate turbine performance losses due to blade LEE.

Sensitivity of BL transition to damage geometry
The aerodynamic performance of rotor blades and, thus, power and AEP of wind turbines hit by LEE strongly depend on the blade BL state [17]. More specifically, the turbine performance depends on whether and to which extent LEE-induced airfoil shape alterations alter the chordwise position of the laminar-to-turbulent BL transition. The dependence of BL transition of eroded blade sections on the freestream Reynolds number and the erosion depth is examined and discussed below, also with the aim of justifying the choice of the fully turbulent flow model for generating the ALPS databases of this study.
The assumption of fully turbulent BLs in eroded airfoil aerodynamics relies on the flow perturbation caused by LEE-induced geometric irregularities being sufficient to trigger BL transition at the LE. However, this is not always the case. At low positive AoA, for example, the flow acceleration on the SS in the LE region is relatively low and the BL is thicker than at higher AoA. In this circumstance, the flow perturbation due to LEE may be insufficient to trigger transition [17].
To highlight the dependence of rough wall BL transition on freestream Reynolds number and roughness height h k , it is convenient to introduce the roughness Reynolds number Re k , given by: where ρ and μ denote fluid density and viscosity, respectively. The variable Re k is sometimes denoted by h + k , the normalized roughness height. The friction velocity u k is defined as: where τ w is the wall shear stress. Introducing the skin friction coefficient , with U ∞ being the freestream velocity, and considering the BL on a rough wall of length c, Eq. (3) can be written as: Experimental evidence indicates that if Re k exceeds a critical threshold (Re k ) crit , that can also be determined using well established correlations [36], transition will occur. The critical roughness height (h k ) crit , above which transition occurs, is thus Eq. (6) shows that (h k ) crit decreases as Re increases. This dependence also highlights that the levels of (h k ) crit at the high Re values at which the outer blade sections of wind turbines operate are relatively small with respect to the airfoil chord c.
To estimate the airfoil BL state for the LE damage type and range considered in this study, the analysis of the aerodynamic performance of the NACA 64-618 subjected to LE damage by moderate size erosion cavities and delamination at Re = 7 M is considered. Each airfoil variant  Fig. 7b. One notes that the transitional and fully turbulent force coefficients of the nominal airfoil differ significantly from each other: the transitional c d is lower at all AoAs, particularly at the higher values of this parameter, and the transitional c l is higher at the higher AoAs, as expected. For given LE damage, conversely, the force coefficients obtained with the transitional and fully turbulent analyses are very close to each other, indicating that at Re = 7 M the flow perturbation due to these damages triggers transition at the LE, resulting in the transitional analysis also predicting fully turbulent BLs on both SS and presure side (PS) for all considered AoA values. With reference to Eq. (6), this occurrence indicates that at Re = 7 M the level of roughness height h k associated with the considered damages is larger than (h k ) crit . This is the reason why fully turbulent simulations have been used for generating all aerodynamic force databases of this study.
The results of Fig. 7 also confirm that the performance loss caused by relatively low erosion levels is larger than that associated with fully turbulent BLs on the nominal airfoil. This confirms the necessity of resolving the eroded LE geometry as done in the generation of the ALPS databases of this study.

Verification of artificial neural network predictions
The construction and verification of the c l and c d ANNs for the erosion damage by LE delamination of the NACA 64-618 airfoil are reported in [20]. The validation of the c l and c d ANNs for the erosion damage by LE cavities are presented below. As stated in Section 3, two separate ANNs for predicting c l and c d of the pitted NACA 64-618 (cavities with Φ/c ≤ 0.002) and the gouged NACA 64-618 (cavities with Φ/c > 0.002) are generated.
The prediction accuracy of the four ANNs generated as explained in Section 3 and based on the databases defined by the geometric parameters of Table 1 is verified by considering a 23-element pitted geometry test set and a 31-element gouged geometry test set, both defined in Table 2.
For each geometry of the two sets, the verification is performed by computing the root-mean-square error (RMSE) of the ML and CFD pre-dictions of the c l and c d coefficient. The average RMSEs of c l and c d for the pits test set are 0.0047 and 2.8210 × 10 − 4 , respectively, whereas the corresponding values for the gouges test set are 0.0242 and 3.7673 × 10 − 4 , respectively.
The outcome of the analysis of the ANNs predictive capabilities for the pitted airfoil is summarized in Fig. 8. In each subplot the CFD and ANNs estimates of the force coefficient of the eroded airfoil are plotted against the AoA, and the CFD estimate of force coefficient of the nominal airfoil is also reported for reference. The two top plots refer to the geometries with the minimum c l and c d RMSE, whereas the two bottom plots refer to the geometries with the maximum c l and c d RMSE. Fig. 9 is structured as Fig. 8 and it refers to the verification of the c l and c d ANNs for gouge test set. These results demonstrate good predictive capabilities of the adopted ANNs and highlight that RMSEs are always smaller than the differences between the CFD estimates of the force coefficients of the eroded and nominal airfoils, an important prerequisite for the wind turbine AEP loss analysis discussed in Section 6.

Results
Here ALPS is used to determine the power curve degradation and the AEP loss of the NREL 5 MW turbine corresponding to three blade LE erosion damages characterized by increasing severity to reflect the progression of the turbine performance loss between consecutive maintenance cycles. The first part of the section provides the definition of the three LE erosion damages; the second and third parts focus on the dependence of AEP losses on both LEE erosion damage type and extent, and wind characteristics of the wind resource, and the fourth part presents a sensitivity analysis to quantify the uncertainty of AEP loss estimates due to uncertain airfoil force data.

Damage definitions
The first LE erosion damage consists of cavities distributed on the  blade LE from 70% rotor radius to the blade tip. The definition of the geometry of this damage is obtained by combining computed data [23][24][25] and available measurements based on laser scans of utility-scale eroded wind turbine blades. These data have been used to build a mathematical function of the eroded volume distribution along the blade span. This distribution is then mapped onto discrete series of pitted and gouged blade strips of the same type as those making up the ALPS database defined in Section 4. The geometric parameters of each strip are set by enforcing consistency with the expected eroded volume at the spanwise position of the considered strip, and varied using random number generators within the definition bounds of the database. The pitted and gouged blade portion extending from 70% rotor radius to the blade tip counts 1930 blade strips. The maximum and minimum values, the mean μ and the standard deviation σ of the five geometric parameters defining the generated LE erosion damage are reported in Table 3.
To define the second LE erosion damage, it is assumed that as erosion progresses, pits and gouges coalesce and initiate LE delamination. In this delamination damage, each of the 1930 strips affected by LE cavities maintains its initial width k 2 Φ and damage extent s l and s u around the LE. The cavities of diameter Φ, however, are replaced by a constantdepth delamination over the curvilinear lengths s l and s u of depth Φ, i. e. one sets d/c = Φ/c. Since the delamination ANNs used in this study are those previously generated and reported in [20], the combinations of the three-parameter sets of delamination parameters obtained as defined above have been suitably filtered to avoid any parameter value falling outside the bounds of the training set.
The third erosion damage, which is the most severe, is also a LE delamination damage, and is defined from the aforementioned milder delamination. The third damage is obtained by the initial damage by erosion cavities adopting the same procedure used to obtain the mild delamination damage but setting d/c = 4Φ/c. The maximum and minimum values, the mean μ and the standard deviation σ of the three geometric parameters defining the moderate (case 1) and severe (case 2) LE delamination damages are reported in Table 4.

Power curve degradation
The ALPS assessment of the performance reduction of the NREL 5 MW turbine featuring the three levels of LEE defined above is performed using OpenFAST and the AERODYN BEMT code [37]. The nominal turbine uses nominal airfoil geometries, with c l and c d data obtained with transitional CFD analyses. In all simulations, the entire turbine structure and its foundations are assumed to be rigid. The assessment of the turbine performance reduction due to each erosion damage is carried out both in the ideal case of uniform freestream wind, often used in preliminary LEE analyses, and the more realistic case of turbulent freestream wind.
The power control strategy of the nominal NREL 5 MW turbine, which is representative of modern speed-and pitch-controlled utilityscale wind turbines is reported in [19,20], and summarized below. The turbine starts producing power at the so-called cut-in speed V cut− in , its power grows in an approximately cubic fashion until the rated wind speed V rated , and is constant from rated to the cut-out speed V cut− out The wind speed interval between V cut− in and V rated (region 2) has constant tip-speed-ratio (TSR) λ, defined by the ratio of blade tip speed and wind speed V. Thus, in region 2 the rotor angular speed varies linearly with V, and the TSR is that maximizing the power coefficient C P . From V rated to V cut− out (region 3) the rotor angular speed is constant, and the blade pitch control keeps the turbine power at its rated value. Before region 2 (region 1.5), and between regions 2 and 3 (region 2.5), a linear relation between torque and rotor speed is adopted. The control settings of all considered damaged turbines use the same overall parameters of the nominal turbines, i.e. rated aerodynamic power of 5.30 MW, rotor speed between 6.9 and 12.1 RPM, and V cut− out = 25 m/s.
In the uniform wind analyses of turbines with eroded blades, the socalled adaptive power control strategy [20] is also applied. The LE geometry alterations due to erosion lead to variations of the optimal TSR λ opt. and the associated optimal power coefficient (C P ) max. . The adaptive control tracks λ opt. of the damaged turbine in region 2, and this enables   some reductions of the power loss due to LEE. The use of the adaptive control would require altering the parameters of the OpenFAST control routine. In this study, as in [20], it has been preferred to implement the adaptive control by wrapping AERODYN with a MATLAB script using gradient-based optimization to track λ opt. , and Newton's method to determine the blade pitch yielding the target rated power above V rated . The curves of power coefficient C P and thrust coefficient C T of the NREL 5 MW turbine operating in uniform wind with nominal blades, and blades affected by each of the three LE erosion damages are plotted in Fig. 10a. All curves are plotted only up to V = 14 m/s, because above this level the power and thrust of all turbines are equal due to the compensation of the blade pitch control [20]. As expected, the C P reduction increases with the severity of the damage, but the reduction is relatively small for the erosion cavity and moderate delamination damages, whereas a substantially larger drop is noted when moving to the severe delamination case. The use of the adaptive control also results in slight alterations of the width of region 2.5 and λ opt. of the four turbines. The values of the wind speed at the start of region 2.5 (V 2.5 ), V rated and λ opt. of the four turbines in uniform wind are reported in Table 5.
The C P and C T curves of the NREL 5 MW turbine operating in turbulent and sheared wind with nominal blades, and blades affected by each of the three LE erosion damages are plotted in Fig. 10b. The selected level of turbulence intensity (TI) 14% is representative of several offshore wind farm sites. Each point of the C P and C T curves is obtained by time-averaging the results of a 10-min OpenFAST simulation of the turbine operating in turbulent wind with the given TI and varying the mean speed from V cut− in to V cut− out . For each mean wind speed and TI, the space-and time-dependent turbulent wind is generated with the NREL TURBSIM code [38]. Wind shear is accounted for by means of a power law profile with shear exponent γ of 1/7, and the OpenFAST default rotor speed and blade pitch control is used in all turbulent wind simulations. Over most of the considered wind speed range, the C P curves of Fig. 10b are lower than those obtained with uniform wind, but the relative distance between the C P curves with turbulent wind and that between the C P curves with uniform wind is different, which indicates that the quantification of AEP losses due to LEE may be different depending on whether a faster uniform wind model or a computationally more costly but more realistic turbulent wind model is adopted. The cross comparison of the turbulent C T profiles of Fig. 10b and their uniform wind counterparts in Fig. 10a also shows that turbulence increases the mean rotor thrust for V <10 m/s. The detrimental impact of turbulence on turbine power and thrust is examined in Fig. 11a, which compares the curves of rotor power P and thrust T for uniform and turbulent wind conditions. The P and T curves for the case of turbulent wind are obtained taking the mean result of OpenFAST simulations considering wind sequences of 10 min, namely the same simulations used to determine the C P and C T curves commented above. The main effect of turbulence is a growing reduction of the power as the rated wind speed of the uniform wind case is approached, an increase of the actual rated wind speed, and a thrust increase both below and above the rated wind speed of the uniform wind case. The variations of the power curve in turbulent wind can be qualitatively explained as follows. For mean wind speeds at or slightly above the rated wind speed in uniform wind, the negative velocity fluctuations make the turbine produce less power than the rated target for part of the 10-min intervals. The mean power is thus smaller than the rated power. For mean wind speeds slightly below the rated wind speed in uniform wind, the positive velocity fluctuations result in the turbine power being capped by the blade pitch control to the rated power. The lower power corresponding to the negative velocity fluctuations is not fully compensated by higher power produced at the speeds resulting from the positive fluctuations. Thus, the mean power is lower than that in uniform wind.
The three remaining subplots of Fig. 11 consider the three LE erosion damages defined above, and each subplot compares the power and thrust curves of a damaged turbine in uniform and turbulent wind to the corresponding curves of the nominal turbine in uniform wind. For each of the damaged turbines the effect of turbulence on the performance curves is comparable to that observed for the nominal turbine. Crosscomparing the turbulent power curves of the last three subplots to that of Fig. 11a points to relatively small differences from cut-into rated wind speed, with the largest differences occurring for the two delamination damages. As discussed below, however, these seemingly small differences can yield significant AEP losses. The sensitivity of the turbulent thrust curve to the considered damages appears to be relatively small.
The OpenFAST turbulent wind simulations of the damaged turbines  used a time step of 0.01 s. The 10-min turbulent wind simulations were run in the serial queue of the HEC cluster, and each of these simulations required an average wall-clock time of 10 h. This runtime is substantially higher than that of conventional OpenFAST analyses because the blades with LE erosion feature about 1930 strips.

Parametric analyses of AEP losses
The impact of using either a uniform or a turbulent wind model in the assessment of AEP losses due to LEE of increasing severity is analyzed in Fig. 12. The figure provides the AEP in MWh of the nominal NREL 5 MW turbine and variants with increasing level of LEE in uniform and turbulent wind. The histograms labeled 'trb. bls.' refer to a turbine with nominal airfoil geometries, but c l and c d data obtained with fully turbulent CFD analyses. This configuration corresponds to the initial stage of LEE, when the LE roughness height reaches the critical roughness height (h k ) crit triggering BL transition at the LE. The successive three damages are those due to LE erosion cavities ('dam. cav.'), and the two LE damages by delamination ('dam. del. 1' and 'dam. del. 2') described above. The white histograms indicate the turbine AEP calculated using the power curve referring to uniform wind, and the black histograms indicate the AEP calculated using the power curve referring to turbulent wind. In all cases, the AEP is determined by integrating the power curve  against an annual wind frequency distribution defined by a Rayleigh distribution with mean wind speed V mean of 9.36 m/s. Except for variant 'trb. bls.', all power curves used for calculating the AEP values of Fig. 12, are those of Fig. 11 with the turbulent wind power curves referring to TI = 14%.
The two expected main trends emerging from the results of Fig. 12 are that a) for a given wind model, the turbine AEP decreases as the severity of the LE erosion damage increases, and b) for given blade geometry the turbine AEP is lower when using the more realistic turbulent wind condition. The former occurrence is due to the fact that, as erosion progresses, c l levels decrease and c d levels increase, reducing the turbine power for V cut− in < V < V rated . For the nominal turbine, the fact that the turbine AEP is lower in turbulent than in uniform wind is caused by the turbulent fluctuations reducing the turbine power around the rated wind speed in uniform wind, due to the power 'filtering' of the blade pitch control, as explained in the previous subsection. For the eroded turbines in turbulent wind, the AEP reduction is aggravated by the poorer aerodynamic perfomance of their eroded blade sections. The percentage AEP difference of the nominal turbine in uniform and turbulent wind is 2.17%, whereas that of the turbine with blades affected by advanced LE delamination (dam. del. 2) in uniform and turbulent wind is 2.49%.
For both the cases of uniform and turbulent wind, the incremental AEP reduction as erosion progresses is particularly pronounced when passing from the LE with erosion cavities to the two configurations with LE delamination. This is because the delamination damages lead to the largest reductions of c l and the largest increases of c d at high AoA. Fig. 12 also shows that, for both uniform and turbulent wind, the AEP of the turbine with pitted and gouged blades is only slightly smaller than the AEP of the turbines with nominal blade geometry and fully turbulent BLs. It is noted, however, that, although small in the considered example, the increase of the AEP loss due to erosion cavities over that of the nominal airfoils with BLs tripped at the LE is significantly affected by the number, position and size of the cavities, and can be larger than in the present case [15]. The dependence of the LEE-induced AEP loss on installation site characteristics is examined in the three subplots of Fig. 13; turbulent wind conditions are considered in all cases. The sensitivity of the AEP loss to V mean is analyzed in Fig. 13a. All sites are characterized by TI = 14% and wind resource defined by a Rayleigh distribution, which can be seen as a Weibull distribution with shape factor k = 2. AEP losses are measured with respect to the AEP of the nominal turbine working at the same sites, and these values are provided in the second column of Table 6 to enable the estimate of the AEP loss of the damaged turbines in MWh. The data of Fig. 13a highlight that, for given LE erosion damage, Fig. 13. Dependence of LEE-induced AEP loss on wind characteristics of installation site: a) sensitivity to mean annual wind speed V mean , b) sensitivity to turbulence intensity TI, and c) sensitivity to shape factor k of Weibull wind frequency distribution.
the percentage AEP loss decreases as V mean increases. This is because increasing V mean the median value of the distribution moves towards the high wind region, where the power reduction due to erosion can be compensated by the blade pitch control [20]. It is noted, however, that even the seemingly modest AEP loss of 0.54% corresponding to damage by erosion cavities of the turbine operating at the site with V mean = 9 m/ s constitute a significant revenue loss. Another important note is that the monetary revenue loss for given damage and varying V mean may have patterns differing from those in Fig. 13a. For example, for the severe delamination damage, the percentage AEP losses for V mean of 7 and 9 m/s are 4.01% and 2.83%, respectively. Using the data of the second column of Table 6, however, the AEP losses associated with the same values of V mean are 567 and 596 MWh, respectively. Thus, although sites with lower V mean incur higher percentage AEP losses, the actual AEP losses differ less.
The sensitivity of the AEP loss to TI is analyzed in Fig. 13b, which considers TI values of 14, 16 and 18%. All sites are characterized by a wind frequency Rayleigh distribution with V mean = 9.36 m/s. AEP losses are measured with respect to the AEP of the nominal turbine working at the same sites, and these values are given in the fourth column of Table 6. The turbine operating at TI of 18 and 16% correspond, respectively, to turbines of class IA + and IA of the International Standard IEC61400 of the International Electrotechnical Commission. The data of Fig. 13b highlight that, for given LE erosion damage, the percentage AEP loss is fairly insensitive to the variations of TI in the considered range. Unlike the case of the AEP loss variation with V mean , the AEP loss diagram using energy units is very similar to that of Fig. 13b, because the AEP of the nominal turbine is fairly insensitive to the considered TI values, as seen in the fourth column of Table 6.
The sensitivity of the AEP loss to the shape factor k of the Weibull annual wind frequency distribution of the installation site is analyzed in Fig. 13c. All sites are characterized by V mean = 9.36 m/s and TI = 14%. AEP losses are measured with respect to the AEP of the nominal turbine working at the same sites, and these values are provided in the sixth column of Table 6. The data of Fig. 13c show that, for given LE erosion damage, the percentage AEP loss increases as k increases. This is because for low values of k the Weibull distribution has high skewness, is biased to low speeds, and its median is lower than its mean. The turbine produces large amounts of power at low speeds, where the percentage AEP loss is relatively low. As k increases, the distribution tends to a symmetric pattern, and median, mean and modal values tend to coalesce. Thus, the turbine produces more power at higher speeds, where the percentage AEP loss is higher. At the highest k values, wind speeds are highly clustered about V mean , with less power produced above V rated and less scope for the pitch control to compensate for the power loss caused by LEE. The sixth column of Table 6 also shows that the AEP of the nominal turbine increases with k. This is both because the considered V mean is fairly close to V rated , and because higher values of k result in a reduction of the time spent at low wind speed. Using these reference AEP data, the AEP losses in MWh associated with severe delamination damage for values of k of 1.5, 2.5 and 3.5 are respectively 432, 715 and 913 MWh. The AEP increase of the nominal turbine increasing k from 2.5 to 3.5 is 1054 MWh. These results suggest that at sites with high values of V mean , high k values are desirable, in general. Taking LEE into account, however, high values of k may result in significant AEP losses, outweighing the benefits of a large k.

AEP sensitivity to uncertainty on airfoil aerodynamics
To estimate the uncertainty level affecting AEP loss estimates due to discrepancies of measured and computed variations of the force coefficients of nominal and delaminated airfoils, the AEP analysis of a variant of the NREL 5 MW turbine is considered. In the modified turbine, the 18%c thick NACA 64-618 airfoil, used by the baseline turbine from 70% rotor radius to blade tip, is replaced by the NACA 63 3 -418 airfoil, which has the same thickness. At a site characterized by TI = 14% and Rayleigh wind frequency distribution with V mean = 9.36 m/s, power curve and AEP of the considered turbine variant are determined using the measured and computed force coefficients of the nominal and delaminated NACA 63 3 -418 airfoil plotted in Fig. 6. The obtained AEPs are reported in the second row of Table 7. The corresponding AEP losses, reported in the third row, show that the differences between measured and computed force data discussed in Section 5.1 result in the airfoil CFD data-based AEP loss estimate exceeding by about 1% of the nominal AEP the estimate based on experimental data. This significant difference underlines the necessity of further investigating and improving the estimate reliability of measured and/or computed aerodynamic force coefficients of delaminated wind turbine airfoils.
To estimate the uncertainty level affecting AEP loss estimates due to discrepancies of measured and computed variations of the force coefficients of nominal and gouged airfoils, the AEP analysis of a second variant of the NREL 5 MW turbine is considered. In this variant, the NACA 64-618 airfoil of the baseline turbine is replaced by the DU 96-W-180 airfoil, which also has thickness of 18%. At the same site considered for the NACA 63 3 -418 turbine variant, power curve and AEP of the second turbine variant are determined using the measured and computed force coefficients of the nominal and gouged DU 96-W-180 airfoil plotted in Fig. 6 of [17]. The obtained AEPs are reported in the fifth row of Table 7. The corresponding AEP losses in the sixth row show that the differences between measured and computed force data observed in [17] results in the airfoil CFD and measurement-based AEP losses differing by only about 0.1% of the nominal AEP. Table 7 indicates that, on the basis of measured airfoil force coefficients, the AEP loss of the gouged DU 96-W-180 turbine variant is greater than that of the delaminated NACA 63 3 -418 turbine variant. This may appear counter-intuitive. This occurrence is partly due to the fact that the wind tunnel measurements of the DU 96-W-180 airfoil are at Re = 1.5 M, whereas those of the NACA 63 3 -418 airfoil are at Re = 5 M. At this higher value, BLs are thinner and less prone to separation, resulting in lower profile losses. The beneficial effect of the higher Re may reduce the performance loss due to LE delamination, expectedly larger than that due to LE gouges for given Re, explaining the observed loss trends.

Conclusions
This study presented an experimentally validated simulation-based Table 6 AEP of nominal turbine for all mean annual wind speed V mean , turbulence intensity TI and shape factor k considered in the analyses of Fig. 13 investigation into the impact of blade LEE on wind turbine AEP. The analyses were based on ALPS, a technology combining the predictive capabilities of RANS CFD with the large runtime reductions enabled by ML. The considered damages by LEE mimicked the progression of erosion in utility-scale wind turbines, and their geometry was generated using available laser scans of eroded blades and previously performed simulations of wind turbine blade LEE yielding estimates of the resulting material loss. The analyses were based on the NREL 5 MW reference turbine, assuming that the outer 30% of its blades were affected by LEE. The mildest damage corresponded to the formation of LE erosion cavities. The next two erosion stages, progressively more advanced, corresponded to the onset of LE delamination and a further degradation of the initial LE delamination damage. Using the turbulent wind model, sensitivity analyses of the AEP of all turbine configurations to the annual mean wind speed, the turbulence intensity and the shape factor of the Weibull wind frequency distributions were performed. Percentage AEP losses varied from between 0.3 and 0.8% of the nominal AEP for the LE damage by erosion cavities to between 2.2 and 4% for the advanced LE delamination. Erosion losses were found to vary fairly little with TI, and significantly more with the mean wind speed and the shape factor. High values of this parameter increase the AEP of the nominal turbine but also result in large AEP losses in the presence of erosion, an occurrence which may outweigh the benefit of a high scale factor.
A comparison of turbine AEPs computed for given nominal or eroded blade geometry, mean speed and shape factor of the wind frequency distribution, using either a uniform or a turbulent wind model found that turbulence reduces turbine AEP by about 2 to 2.5%. A sensitivity analysis of the dependence of the LEE-induced AEP loss on variations of the airfoil force coefficients pointed to low uncertainty of about 0.1% of the nominal AEP for AEP losses due to erosion cavities, whereas a higher uncertainty of about 1% of the nominal AEP was found for the case of LE delamination. This indicates the necessity of further investigating the accuracy of CFD and experimental methods for delaminated airfoil aerodynamics and reducing the uncertainty on the force coefficients of these airfoils.
A sensitivity study of the impact of LEE on BL transition and aerodynamic performance showed that real LE erosion features, which are sparser than the distributed roughness often considered in simulations and experimental tests, and have mean roughness height larger than that needed to trip BLs at the LE, induce aerodynamic losses well above those caused by the loss of BL laminarity. For this reason, it is advisable to geometrically resolve these features rather than modeling their effects in eroded blade aerodynamics.
The generalization of ALPS to LE erosion cavities required automating geometry and mesh generation, and CFD analysis and postprocessing of a large number of 3D eroded blade strips, whose force coefficients are required for building the databases to train and validate the adopted ANNs. The ALPS use of BEMT codes enables very high radial resolution of the LE damage geometry (nearly 2,000 blade strips were used over the outermost part of the NREL 5 MW turbine blades in this study), but further work is required to enable the analysis of LE erosion phases characterized by hybrid damage patterns, such as non-circular grooves and delamination patterns with ragged contours.
The CFD-based generation of the ALPS databases for damage by LE cavities and delamination requires about two weeks of medium-size high-performance computing. Thereafter, ALPS can determine the power of hundreds of turbines whose blades feature LEE in just a few minutes. Thus, the technology is a key tool to support wind farm cost analysis and, ultimately, implement predictive maintenance strategies.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.