Impact of meteorological data factors and material characterization method on the predictions of leading edge erosion of wind turbine blades

Leading edge erosion of wind turbine blades is a major contributor to wind farm energy yield losses and maintenance costs. Presented is a multidisciplinary framework for predicting rain erosion lifetimes of wind turbine blades. Key aim is assessing the sensitivity of lifetime predictions to: modeling aspects (material erosion model, blade aerodynamics), input data and/or their preprocessing (joint frequency distribution of wind speed and droplet size based on synchronous site-specific measurements versus frequency distribution generated with partly site-agnostic modeling standards, wind speed records of nacelle anemometer or extrapolated at hub height from met masts), and environmental conditions (UV weathering). The analyses consider a Northwest England onshore site where a utility-scale turbine is operational, focus on a reference 5 MW turbine assumed operational at the site, and use a typical leading edge coating material. It is found that the largest variations in erosion lifetime predictions are due to material erosion model (based on rain erosion test data or fundamental material properties) and wind and rain model (measurement-based joint wind speed and droplet size distribution or standard-based modeled distribution). The use of joint wind and rain distribution also enables identifying wind/rain states with highest erosion potential, knowledge paramount to deploying erosion-safe turbine control.


Introduction
Rain erosion on wind turbine (WT) blades spoils the aerodynamic performance of the rotor, leading to energy yield losses and increased maintenance costs of wind farms.Blade leading edge erosion (LEE) may result in annual energy production (AEP) loss ranging from 1% to 4% [1][2][3].The progression of LEE is influenced by several factors, including the coating's resistance to fatigue caused by repeated impact with rain droplets, weathering, and the conditions of the impinging rain droplets, characterized by number of impacts, droplet size and impact velocity.Therefore, a rain event can lead to varying LEE degrees, depending on wind speed, blade size and material, and rotor speed.For instance, heavy rain at low wind speeds yields no erosion, since the rotor is idle and the impact speed of the droplet is fairly low.Identifying wind and rain conditions critical to LEE is necessary to define sitedependent strategies for mitigating LEE, and requires reliable LEE assessment capabilities.
Reliable LEE predictions require combining models to account for all factors affecting the phenomenon, yielding a multiphysics and multidisciplinary methodology subjected to various sources of uncertainty.Some of these uncertainty factors, such as those associated with the droplet size distribution (DSD), have been recently examined [4][5][6].These studies quantified the sensitivity of rain LEE assessments to disdrometer type, DSD and rain flux modeling, and highlighted how standardized parametric DSDs, like the Best [7] and the Marshall and Palmer [8] DSDs can introduce significant errors in LEE predictions.It is essential, however, not only to use DSD sufficiently representative of the site under consideration, but also to provide the rain LEE prediction system with concurrent wind and rain data, i.e. a joint frequency distribution function (FDF) of wind speed and droplet size.This is one of the problems addressed in the present study, which uses a joint FDF of wind and rain data based on synchronous anemometer and https://doi.org/10.1016/j.renene.2024.120549Received 2 February 2024; Received in revised form 18 April 2024; Accepted 22 April 2024 disdrometer measurements carried out at an onshore site in Northwest England over one year.
The meteorological characterization of wind farm sites is one of the key aspects of a LEE prediction framework.Recent studies presented comprehensive approaches for LEE damage assessment, an essential input for developing optimal site-specific LEE mitigation strategies.Verma et al. [9,10] used measured rain data to define a statistical rain model which, along with measurement-based wind speed probability distribution functions (PDFs) and Springer erosion model [11], enabled finding distributions of critical rain parameters for LEE assessment.The model was applied to assess the LEE potential of different sites in the Netherlands, highlighting the sensitivity of the coating material life prediction to atmospheric conditions.Letson et al. [5] adopted an energy-based erosion model to estimate the damage accumulation for a variety of DSDs.Prieto et al. [12] proposed an analysis framework that integrated meteorological data, operational settings, and Springer model to characterize LEE severity, taking into account the effect of rain, snow, sea spray, and fog.This framework was validated using data from five wind farms, and a good qualitative alignment with observed LEE patterns was obtained.Visbech et al. [13] employed a data-driven framework for modeling erosion damage, based on blade inspections from several wind farms in northern Europe and a mesoscale numerical weather prediction model, thus deriving erosion maps over the considered region.
In most of the above studies, erosion damage modeling was centered on the occurrence of coating material failure at the leading edge (LE), utilizing corrections to account for the effect of blade geometry, or simply assuming a normal impact to ascertain the time to the damage onset on a generic point of the LE.
Recently, researchers have proposed approaches to more accurately resolve the characteristics of the impact of rain droplets on the blade by computing the distribution, velocity and angle of impact considering the local blade geometry and the effect of aerodynamics on rain droplet impact trajectories.Herring et al. [14] used Computational Fluid Dynamics (CFD) and particle tracking simulations to model the impact characteristics of a droplet, the number of impacts, and the effect of aerodynamics in swirling arm rain erosion test rigs.Using swirling arm rain erosion testing (RET) data, the study highlighted that aerodynamic effects are significant when the droplet size exceeds certain thresholds.CFD and particle tracking methods for modeling the impact patterns on WT blades have also been widely used by Castorrini et al. [15,16].
Verma et al. [17] demonstrated how other site-dependent factors, such as wind turbulence and rain intensity, influence the rate of LEE.Other studies focused on droplet/material interaction modeling, and the erosion model itself [18][19][20].A review of the main issues and approaches is reported in [21].
Consideration of the actual blade geometry and the aerodynamic effects on the impinging droplet trajectories in determining the droplet impact characteristics enhances the accuracy of LEE predictions.This is because the impact angle affects the normal component of the impact velocity, and the impact angle is defined by the relative velocity of the droplet and the normal to the blade surface at the impact point.A common approach to account for aerodynamic effects is defining the so-called impingement efficiency, an empirical correction of the impact number based on simulations and experiments [22,23].The studies of [15,16,24] computed the impact pattern using high-fidelity CFD and particle tracking simulations to account for blade geometry and aerodynamic effects.In [25] the high-fidelity simulation approach was used for blade LEE analyses in the general case of variable wind, rain and WT operating conditions.The cost of these comprehensive analyses was reduced by creating a database to support the generation of computationally efficient surrogate models of the impact pattern.
One fundamental component of the LEE analysis framework is the calculation of the rain erosion damage of the blade coating material.Models for this have been developed primarily from Springer model [11], an approach based on fundamental material properties obtained from fatigue and ultimate strength tests.When applied to WT blade coatings, however, this approach was found insufficiently accurate [26].One of the intermediate steps of LEE analyses requires using the so-called V-N curves.A V-N curve relates the droplet impact velocity to the number of impacts at which the erosion damage onset is detected.This curve can be determined from the fundamental material properties, as in Springer model, or measured with RET.However, the analyses of [26], considering different LE coating materials, highlighted notable discrepancies between the V-N curves measured with RET and those computed using the fundamental properties.
Several modifications and alternatives to Springer model have been proposed.Notable examples include the Siemens model [27], the TNO model [18,28], based on fatigue rules, and the DTU model [29] based on the kinetic energy of the impinging droplet.Moving from the approach of [27], DNV published the DNV-GL-RP-0573 report [30] outlining a standardized framework for performing rain erosion analyses based on RET data.
A swirling arm rain erosion test is typically accelerated and performed at single impact velocity and droplet size.In field operation, however, there is significant variability of droplet size and impact velocity, the former parameter being site-specific and the latter parameter depending mainly on blade geometry and WT operating condition.Various approaches have been proposed to account for this variability.Some studies handled the rain variability quantifiable by DSDs and rain flux, applying corrections based on droplet size or impact energy [5,9,26] Bech et al. [31] presented swirling arm RET results for different droplet sizes, highlighting the limitations of both the approaches based on the baseline Springer model and the energy-based model.The study also proposed a droplet size-dependent empirical model yielding the erosion damage as a function of the impact velocity.
The work herein critically analyzes and brings together most of the techniques outlined above into an integrated blade erosion prediction system.This framework is used to examine the impact that different factors and model choices have on LEE assessments.The framework combines modules for considering wind and rain distributions, WT operating conditions, droplet impact patterns, damage incubation and propagation.Using the field measurement-based joint DSD and wind speeds, the framework also enables reliably predicting the coating lifetime, and identifying the combined wind and rain regimes yielding maximum erosion rate.Presented parametric analyses assess the impact of different model choices on LEE damage predictions.Considered modeling aspects are: (i) inclusion of aerodynamics in the computation of droplets impact patterns, (ii) adoption a method for damage and lifetime assessment based on either RET data or fundamental material properties, (iii) inclusion of weathering effects in LEE analyses, and (iv) source of rain and wind data, i.e. use of data based on either recommended DSD and wind models or measurement based wind speed and droplet size joint FDFs.
The article is organized as follows.The LEE analysis framework is presented in Section 2. The model set-up and the considered case study are defined in Section 3, made up of four subsections.Section 3.1 summarizes the two main approaches for rain erosion modeling, namely that based on fundamental material properties, i.e.Springer model, and that based on RET-inferred material properties [26,30].Furthermore, material properties from RET results are derived for the considered coating material with and without weathering.Section 3.2 presents the derivation of the wind speed and droplet size joint FDF at the site considered in this study, whereas Section 3.3 describes the considered WT and the blade geometry representation used for the analyses.Section 3.4 presents the modeling approach for computing the droplet impact patterns with and without aerodynamic effects.The results are presented in Section 4. These are based on the erosion analysis of the considered WT assumed to operate at the aforementioned Northwest England onshore site, and include comparative analyses of damage incubation time and progression.The main findings of the investigation are presented in Section 5.

Blade LEE prediction system
The multidisciplinary blade LEE prediction system incorporates several modules.Its workflow is presented in Fig. 1.The analysis starts with four pre-processing operations, namely (1) definition of material properties and relevant quantities required by the erosion model, (2) set-up of atmospheric, i.e. rain and wind, data, (3) blade geometry processing and (4) generation and organization of blade aerodynamic.
Material properties are determined in accordance with the specifications and requirements of the rain erosion model.In this study, results obtained with two alternative approaches to erosion modeling are compared.One method is the original Springer model [11], based on the fundamental properties of the coating and substrate materials determined with conventional structural tests; the other method is the erosion model of the DNV-GL-RP-0573 standard [30], based on RET data.Both approaches are explained in detail in Section 3.1.The investigations of this study also use the erosion model of [30] to compare the erosion characteristics of a LE coating material with and without UV protection to assess the impact of weathering.
Rain and wind data are processed to generate the FDF   (  , ), where   is the wind speed and  is the droplet size in mm.The FDF enables obtaining the number of droplets per unit surface with diameter in the bin ( − , ) impacting the blade for wind speed is in the bin (  −   ,   ).Typically, the yearly rain and wind data distributions for erosion analyses are generated with independent analytical representations of wind velocity, rain flux and droplet size distribution [30].In the present study, erosion failure results obtained using the standard approach to wind and rain modeling [30] are compared to those obtained using a joint FDF   (  , ) based on anemometer and disdrometer data.The procedure for rain and wind data acquisition at the considered onshore site, and their processing are described in Section 3.2, which also describes the generation of   (  , ) based on the DNV standard [30].
The computation of the blade damage uses the approach of [32], which is in line with the damage and durability definitions of [30], but considers also the blade geometry and local aerodynamics.The blade is subdivided into chordwise strips, and the erosion damage is calculated at all points of the mean profile of each strip of mean chord .The computed erosion damage depends on the normalized curvilinear coordinate ∕, which is 0 at the LE, negative on the suction side and positive on the pressure side of the strip.For each blade strip and mean wind speed, a dataset of relative turbulent aerodynamic states, which is also affected by WT control dynamics, is computed.The key operating variables are the strip relative angle of attack (AoA)  and flow velocity   .For each   , the strip operating variables are grouped in a table with entry {(,   )  }, with  = 1, … ,   , which collects all the   records.The selected WT and its analysis code, the blade strip set-up and the description of the numerical method to build the {(,   )} tables are described in Section 3.3.
Following the pre-processing phase described above, the incubation time of the erosion damage is calculated by executing the five nested loops indicated with ellipses in Fig. 1.The first loop iterates over the blade strips.The second loop runs over the   {(,   )  } entries for each mean wind speed   .The third and fourth loops, merged in a single process in the workflow of Fig. 1, iterate over the   mean wind speed and   particle diameter bins of the joint FDF   (  , ).The final fifth loop scans all the all   segments defining the strip airfoil.For a given strip, the average damage  in year −1 of the th segment is computed as in [30] by means of a Palmgren-Miner rule summation giving: where indices  and  refer, respectively, to particular values of droplet diameter   and impact velocity   , and index  indicates the th airfoil segment of the blade strip under analysis.For each triplet (, , ),   is the number of droplets per unit surface impacting the blade that would give  = 1, whereas   is the number of droplets per unit surface that actually impact the blade.The value of   depends also on the material and erosion model, discussed in Section 3.1, whereas that of   depends also on the wind and rain FDF   (  , ), the blade geometry model and the aerodynamic effects, when included.Its calculation is described in Section 3.4.
The coating material durability is thus given by: After the damage incubation threshold (   = 1) is reached at the th segment of a blade strip, the impacting droplets begin to remove material from that segment.Following the approach in [16], the erosion depth  normal to the ℎ segment is calculated incrementally by accumulating the term: for each aerodynamic state (,   )  following the end of incubation, where   is the mass per unit surface removed by the   droplets and   is the coating material density in kg/m 3 .The term  is calculated by looping over all droplet sizes and mean wind speeds, as done for the damage  defined by Eq. (1).To calculate   , the approach presented in [16] is adopted.The method, which builds upon Springer's original work [11], gives: where   is the droplet mass in kg and   is the droplet diameter in m.
In this work, it is assumed that erosion stops when  reaches the total coating thickness.Erosion of the substrate materials is not considered.

Definition of models set-up and test cases
This section defines the models used for all the pre-processing stages summarized in the previous section.A test case is also defined for each modeling set-up.The results obtained using different model choices in terms of damage incubation  or, equivalently, durability 1∕, and erosion depth and pattern  will be reported and compared in Section 4.

Erosion models
The LE protection (LEP) system considered in all analyses of this study is a prototype polyurethane coating named ProtectedUV_ISO 16474-2_0h_V90_Scotland, developed by Aerox Advanced Polymers.Part of the erosion characteristics of this LE coating have been derived from swirling arm RET data performed in accordance with the DNV-GL-RP-0171 standard [33].The aforementioned coating was tested in the swirling arm rain erosion tester of AeroNordic and Polytech.The rigs are based on the wind industry design guidelines of the DNV-GL-RP-0171 standard for testing LEP systems.Table 1 reports the parameters of the considered tests.Three different measurement sets were performed, each with a different angular speed of the swirling arms.Sets 1 A and 1B, using the AeroNordic facility, refer to the coating without weathering, whereas set 2, using the Polytech dacility, refers to the weathered coating, as described in Section 3.1.3.The structure and geometry of the test specimens are shown in Fig. 2. Each specimen has spanwise length of 39 mm and its center is at 1 m from the rotational axis.
The fundamental material properties and fatigue properties of the coating are given in Table 2, whereas the properties of the materials making up the LEP system are reported in Table 3, Fig. 3 shows one specimen of test set 1 A at the start and the end of the incubation time,   whereas the two images of Fig. 4 depict the initial and final states of a weathered specimen of test set 2.
As mentioned above, this study also investigates the effect of weathering on the erosion performance of LEP systems.To investigate this effect, the considered LEP material was tested also after exposure to UV radiation (test set 2).In an accelerated weathering treatment implemented in accordance with the ISO 16474-2 standard [34], the samples were exposed for 1000 h to xenon-arc light in the presence of moisture to reproduce the weathering effects.Fig. 5 shows the weathering equipment and the RET specimens.Fig. 6 reports the RET data of all three test sets, i.e. for the case in which the specimens underwent RET without weathering and that in which the samples were tested in the swirling arm facility after 1000 h of weathering.The following subsections describe the two approaches used in this study to determine the erosion characteristics and lifetime of material coatings, with one of the two approaches requiring the type of data reported in Fig. 6, and also define four frameworks for two comparative case studies aiming to establish the impact of the erosion lifetime assessment and weathering on the coating lifetime.

Derivation of erosion strength from fundamental material properties
In this section, the method for calculating the coating erosion life making use only of fundamental material properties is presented.This method, due to Springer [11], enables determining the coating durability making use of data like those provided in Tables 2 and 3 for coating ProtectedUV_ISO 16474-2_0h_V90_Scotland. Springer model is reviewed in [26] and summarized below.
During the incubation period, the total number of impacts per unit area  is related to the total number of impacts per site  * (i.e. at a given point of the impacted surface) by where  is the droplet diameter in mm.Springer assumed that the total number of impacts per site in the incubation period can be evaluated with the equation where   is the erosion strength of the coating,  0 is the average stress on the surface of the coating, and  1 and  2 are empirical constant parameters.Springer provided values of  1 and  2 of 7 × 10 −6 and 5.7, respectively.These values were obtained by fitting Eq. ( 5) to the experimental data available at the time the model was developed.
Inserting the values of the constants  1 and  2 and Eq. ( 6) into Eq.( 5) gives: The average stress of the coating surface is defined as where  is the water hammer pressure,   is the relative acoustic impedance of liquid droplet and coating,   is the relative acoustic impedance of coating and substrate, and  is a thickness parameter that accounts for the attenuation of the reflections at the coating layer.The aforementioned impedances are defined as: where  =  is the impedance of the material (subscripts L, C and S refer to the liquid, coating, and substrate, respectively,  denotes density, and  is the elastic wave speed.The thickness parameter may be calculated as    where ℎ  is the coating thickness. If the value of the relative impedance of substrate/coating interface   equals zero, then the coating material is considered to behave as the substrate, and Eq. ( 8) reduces to  0 =  , which may be used to compute the average stress for homogeneous materials.The water hammer pressure can be defined as where  and  are the droplet impact velocity and angle, respectively.Inserting Eq. ( 12) into Eq.( 8), the average stress of the coating surface becomes The erosion strength of the coating   is defined as a function of its material properties, i.e.
where   is the ultimate tensile strength of the coating,   is the endurance limit,   is a Wöhler fatigue constant, and   is Poisson's ratio.The erosion strength   is a suitable comprehensive parameter describing LEP capabilities to withstand the impact energy in fatigue analysis.The final expression of   is justified by the fact that ( 1 and   has typically order of magnitude 10.The value of the ultimate strength obtained in conditions of static testing was found to give overly unsatisfactory results of the erosion characteristics with respect to experimental observations.Therefore, an alternative approach, following the approach of [35,36], is adopted to determine   for the polyurethane coating considered in this study.The dependence of   on the strain rate ε is modeled by means of power law of the form: where  is a material parameter that depends on temperature, strain, and material microstructure, and  is the strain rate sensitivity index.Using a strain rate of 800 Hz in Eq. ( 15), one finds that the values of the two parameters for the considered coating are  ≈ 0.03 and  ≈ 1.01.Further detail on the derivation of the numerical values of these two parameters can be found in [26].
An important issue in fatigue analysis is how to determine the slope   of the fatigue curve, a parameter difficult to obtain experimentally for typical elastomeric materials used for LE coatings.In the range of stress   in which fatigue occurs, the Wöhler fatigue curve, relating stress and material life N, reads: Once the data of the whole Wöhler curve have been determined experimentally, including the endurance limit (  ), the exponent   can be estimated as where and   is the number of fatigue cycles at the knee of the Wöhler curve, i.e. the value above which the curve flattens out.
To evaluate the average stress value at the coating/liquid and coating/substrate interfaces during the impact, the average number of mechanical wave reflections in the coating material , is introduced: The parameter  depends both on the coating thickness ℎ  and the droplet diameter , due to the dependence of  on these two variables, as shown by Eq. (11).Eq. ( 14) is finally restated as an equivalent erosion resistance parameter of the coating, namely: Springer model assumes that if there is no stress wave reflection in the time of duration of the impact pressure at the liquid/coating interface, then the coating material can be assumed to behave as if it had infinite thickness.Springer showed that a coating material can be considered of infinite thickness, i.e. the LEP system treated as a homogeneous material, if: Because the thickness used for LE coatings is small, typically a few tenths of a millimeter, these coatings are typically considered as thin, regardless of the drop size.
In the analyses of this study, four cases are considered for determining the erosion characteristics of the considered coating.Case 1 is that in which   is calculated with Eq. (19), which depends on all the fundamental properties of the LEP system.

Derivation of erosion strength from measured V-N data and performance analysis based on linear regression
This section describes an alternative method for determining the erosion characteristics of the coating, partly based on RET data.The fatigue curve of droplet impact velocity versus number of impacts, known as V-N curve, is derived from RET data.The number of cycles  to end of incubation is the number of droplet impacts, and the stress variable on which it depends corresponds the impact velocity  .In the context of LEE, the V-N curve is described by a power law of the type which, in a log-log graph, becomes a linear law with slope − and intercept log(): Following the ASTM E739-10:2015 guidelines, the experimental data are fitted with a linear regression in a log-log representation.Statistical methods are used to obtain the V-N curve with a confidence level of 95% [37].Using the methodology above, it is found that  = 1.3007 × 10 13 and  = 4.97.Following the guidelines of the DNV-GL-RP-0573 standard [30], once the fatigue slope  is derived, the erosion strength   fit of the LEP system is obtained from Eq. ( 7) by replacing the exponent 5.7 with the fitted exponent : 8.9 where   fit is the number of impact to incubation observed in the RET experiment and The symbol    denotes the impact velocity in the RET experiment corresponding to   fit .
The equivalent raw strength of the coating   [30], which is one of the input parameters of Springer model, is calculated using the mean droplet diameter used in the test: where  and  are obtained for the RET data using Eqs.( 11) and ( 18), respectively.
The parameter   is then used to calculate the values of   for each droplet diameter considered in the erosion analysis: The allowable impacts on the surface per unit area to the end of incubation for any combination of drop diameter and impact velocity is calculated using the above equations for  0 and the value of   given by Eq. ( 26): where the diameter  is in mm.In Section 4, results referring to calculation of the erosion characteristics and coating life based on the method described in this section, are labeled Case 2.

UV-degraded material
As mentioned, this study considers also the impact of UV weathering on the erosion characteristics and lifetime of LE coatings.The methods for determining the erosion characteristics of weathered coatings are the same as those for materials not exposed to weathering, described in Sections 3.1.1and 3.1.2.Predicted lifetimes of coatings without and with weathering differ only because the fundamental material properties and/or the RET-based V-N curves of the two materials differ.In Section 4, the comparative erosion analyses of the coating subjected to weathering by UV exposure and that without weathering are labeled Case 3 and Case 4, respectively.
The intention was to perform the analyses of Cases 3 and 4 using RET based V-N curves, using the method discussed in Section 3.1.2.However, since there were not enough V-N data for the weathered samples, the exponent  for Cases 3 and 4 were assumed equal to the exponent  2 = 5.7 of the baseline Springer model.The RET data for both materials are then used to fit the erosion strength   , as discussed in Section 3.1.2.

Summary
Table 4 reports the values of the erosion strength parameters and the exponent  for the Cases 1-4 defined above.The droplet diameter  of 2.61 mm of the RET experiments is used in all cases.As discussed above, the erosion strength parameters in Case 1 are determined by using Springer model with the fundamental material properties, whereas in Case 2 these parameters are determined by means of RET-data using the DNV-GL-RP-0573 framework.In Cases 3 and 4 the DNV-GL-RP-0573 is also adopted except for the fact that the exponent  is set to the Springer value of 5.7, as discussed in Section 3.1.3.Fig. 7 reports all RET data used in Cases 2-4, along with the V-N curves used for all four cases.The V-N curve of Case 1 was generated using the original Springer model, i.e. making use only of the fundamental material properties.The V-N curve of Case 2 was generated by fitting the corresponding RET data set, whereas the V-N curves of Cases 3 and 4 were obtained by using the theoretical exponent of 5.7 and fitting the available RET data.It is noted that the differences among the fitted V-N curves of Cases 2-4 are significantly smaller than the differences between these three curves and the curve of Case 1 based on Springer model.More specifically, with reference to the RET conditions, Springer model predicts significantly shorter erosion lives than any of the RET-data based curves.The quantitative impact of these differences in the context of real WT operation is analyzed and discussed in Section 4.

Rain and wind data
Synchronous wind and rain data measurements were performed at the Hazelrigg UK Met Office weather station (climatological station 7236) at Lancaster University [38], shown in Fig. 8.The station is at a latitude of 54.0138 • and a longitude of 2.7749 • (grid reference 490 579), 1 km north-east of Lancaster University campus and 10 km east of the Irish Sea coast.The site is at 95 m above sea level, and is surrounded by trees and agricultural land.A 2.3 MW WT is operational about 200 m south-west of the station.The prevailing wind direction is south-west.Data recorded at Hazelrigg have been widely adopted over the last two decades to evaluate the quality of the air in Northwest England in terms of pollutant, dust and particles transport [39][40][41].
At the station, a tipping bucket (TB) rain gauge (model ARG100, Campbell Scientific) is used for rainfall measurements.It provides rainfall data in mm every 10 min.The gauge sensitivity is 0.20 mm of rain per tip, and the maximum rainfall rate it can measure is 500 mm/h.Data on air and ground temperature, air pressure, total and diffuse solar irradiance, drizzling and hail are also available at the station.A cup anemometer (model A100LK-L, Campbell Scientific), enabling wind speed measurements according to the IEC 61400-12 Class 1 standard, measures the wind speed in 10-min averages.The anemometer has a measurement range of 0 to 77 m/s, an accuracy of 1% of the measured speed and is positioned at the top of a 10-m met mast.
The Supervisory Control and Data Acquisition (SCADA) data of the 2.3 MW WT are also available.The WT nacelle is at 90 m above the ground.The SCADA data are collected as 10-min averages.Minimum and maximum values of all measured parameters, including wind speed measured by a cup anemometer on the nacelle, rotor speed, active and reactive power and nacelle yaw angle are also recorded.
From February 2017 to October 2019, a disdrometer of the Disdrometer Verification Network (DiVeN) project [42] was also operational at the Hazelrigg station.The instrument is a Thies™ laser precipitation monitor, which uses an infrared beam with a wavelength of 785 nm, and has horizontal area of capture of 45.6 cm 2 .It measures drop diameters in 20 bins from 0.125 to 8 mm, and droplet speeds in 22 bins from 0.0 to 20.0 m/s.Data are recorded at a rate of 1 min.
The generation of the reference joint FDF   (,   ) of wind speed and rain droplet size used in this study is based on synchronous 356 consecutive-day measurements of wind speed by the WT nacelle anemometer and rain droplet size by the DiVeN disdrometer from October 2018 to September 2019.For only 5 days of September 2019, the disdrometer did not record any data, and this gap in the 356-day reference period, corresponding to 1.4% of the period, was filled in with data referring to the same days of September 2017.
The reference joint FDF   (,   ) distributes the droplet diameter data in 20 bins from 0.125 to 8 mm, and the wind data in 30 bins from 0 to 30 m/s.To synchronize the DiVeN rain data and the wind data of the WT anemometer, an aggregation step of the rain data was required, since the wind data refer to 10-min intervals and the rain droplet data to 1-min intervals.The aggregation consisted of redistributing, i.e. recounting, the rain droplets of each set of 10 consecutive 1-min intervals in the same aforementioned 20 diameter bins.This step was repeated for all groups of chronologically ordered sets of 10 1-min intervals from the start to the end of the 356-day reference period.Following this operation, wind and rain data, whose measurement started at exactly the same time, were both referred to 10-min intervals, and could thus be merged to generate the desired FDF   (,   ).
In order to confirm the reliability of the DiVeN droplet size data collected over the reference period, a cross-comparison of the rainfall rate frequency distributions generated by using the TB data and by integrating the DiVeN droplet data was carried out.In both cases, rainfall rates are computed over 1-h intervals.In the DiVeN case, the rainfall in each 1-h interval (rainfall rate) is computed by adding up the volume of all detected droplets during 60 consecutive 1-min intervals, and dividing the total water volume by the area of the capture section of the disdrometer.In the TB case, the hourly rainfall is simply the sum of the rainfall during 6 consecutive 10-min intervals.The two rainfall rate distributions are compared in Fig. 9, which shows excellent agreement between the two estimates, thus confirming the suitability of the DiVeN droplet data distributions for the blade erosion analysis of this study.The rainfall in the 356-day reference period computed by using the TB-based and the DiVeN-based rainfall rate distributions are 1197.8mm and 1184.3 mm, respectively, which differ by less than 1%.The DiVeN disdrometer also calculates the rain flux for each 1-min interval.A further verification step was performed by calculating the rain fluxes for each 1-min interval using the droplet volumes, similarly to the process above for the 10-min intervals.This result matches the value reported by the disdrometer in about 80% of the 356-day reference period.
As mentioned above, the wind speed at the Hazelrigg site is also measured by an anemometer mounted at 10 m on a met mast.To compare the nacelle and met mast wind records, the latter set of data was extrapolated to the WT nacelle height by means of the wind shear power law: )  (28) using  = 0.2, in accordance with the guidelines of the IEC 61400-1 standard for onshore applications.Then, both wind data records were fitted with a Weibull PDF.The shape factor, scale factor and mean wind speed of the Weibull PDF based on the extrapolated wind data are 1.73, 6.13 m/s, and 5.46 m/s, respectively; the corresponding values for the nacelle-based data are 1.9, 6.8 m/s and 6.03 m/s, respectively.Fig. 10 compares the two Weibull PDFs.In line with the different PDF parameters reported above, indicating a larger mean speed for the nacelle data, the two curves differ by non-negligible amounts.This occurrence may be due to both atmospheric physics effects not properly accounted for by the simplistic representation of the wind shear layer of Eq. ( 28) (e.g.time-dependent air/ground heat transfer unaccounted for by a time-and height-independent exponent  in Eq. ( 28), and complex WT aerodynamic effects (e.g.highly turbulent flow at the nacelle anemometer due to the upstream rotating blades).The differences of the two PDFs in Fig. 10 are significant with regard to AEP assessments.Part of the analyses below aims at verifying if these differences impact significantly also the outcome of LEE analyses.The left map of Fig. 11 is the reference FDF   (,   ), whereas the right map is a modified FDF that differs from the former one because it uses the wind data from the met mast anemometer.The two FDFs are quite close to each other, with some differences visible mostly for wind speed above 20 m/s and droplet size between 1.5 and 2 mm.Both FDFs are used in the analyses below to assess the impact of uncertain wind data on the LEE analysis outcome when using joint wind and rain data distributions.
The DNV-GL-RP-0573 standard provides a methodology for defining the rain load for LEE assessments for cases in which no detailed measured rain data are available.The approach combines statistical models of DSD and rainfall rate.More specifically, the model proposed by Best [7] is used to obtain the yearly PDF of the fractional volume occupied by droplets of diameter  (in mm) for a particular rainfall rate  (in mm/h).This PDF reads: )  (29) where  = 2.25,  = 1.3 and  = 0.232.Eq. ( 29) is then used to obtain the yearly DSD for a given rainfall rate : where  = 67,  = 0.846 and   is the volume in mm 3 of the rain droplet of diameter , assumed to be spherical.It is noted that   is the total yearly amount of water corresponding to rainfall rate , and the PDF   () has dimensions mm −1 .The PDF enables calculating the number of droplets per cubic meter of air which have diameter in a selected interval.
The annual frequency of the rainfall rate I is given by: where  denotes the total annual rainfall in mm, and recommended values for   and   are −0.8 and 1.2, respectively.
Fig. 12 compares the total annual DSD computed with the DNV-GL-RP-0573 approach and that derived from the summation of the DiVen data.The two distribution are in good agreement for droplet size from 0.5 to 2.5 mm.Conversely, the DiVen data-based DSD shows a smaller number of droplets below 0.5 mm, and larger number of droplets with diameter above 2.5 mm.Fig. 13 compares three FDF maps   (,   ) of droplet size and wind speed.The map of Fig. 13(a), labeled 'DiVeN' is the reference joint FDF, and it is the same map reported in Fig. 11(a).The map of Fig. 13(b), labeled 'DNV (Best)', is obtained by using the DNV-GL-RP-0573 approach.The underlying Weibull distribution is that fitting the wind speeds recorded by the WT SCADA system; the total yearly rainfall used in Eq. ( 31) equals that obtained from the DiVen data.For each wind speed   considered in the map of Fig. 13(b), the Best DSD of Fig. 12 is used to calculate the number of droplets of all considered diameters .Thus, the number of hours for which droplets with diameter  are observed is proportional to the number of hours for which the wind speed has value   .The map of Fig. 13(c), labeled 'DNV (DiVeN)', differs from that of Fig. 13(b) only because the former uses the DiVeN annual DSD reported in Fig. 12 rather than the notional Best DSD depicted in the same figure.All three FDFs of Fig. 13 are used in the comparative analyses of Section 4.3 to assess the sensitivity of the predicted coating life to the adopted wind and rain model.More specifically, the FDFs of Figs.13(a) and 13(b) are used to assess the coating life sensitivity to using synchronous co-located wind and rain measurements or the notional model of the DNV-GL-RP-0573 standard.The analysis based om the FDF of Fig. 13(c) aims, instead, at assessing separately the effect of using a notional rather than a measured DSD, and the effect of neglecting the measurement-based wind and rain time-correlation.

Wind turbine model
The analyses below use the NREL reference 5 MW WT [43], which has hub height of 90 m, as the Lancaster University WT described in Section 3.2.The turbine model has been implemented in Open-FAST [44], an aero-servo-elastic and multibody simulation framework for WT simulation.OpenFAST is used to simulate the WT operation from cut-in to cut-out speeds, and extract the time-series of blade strip aerodynamics (time-dependent wind relative velocity and AoA) required for calculating the erosion damage.The WT control uses variable rotor speed from cut-in to rated wind speed of 11.5 m/s to maximize power, and variable blade pitch to maintain power at the rated level from rated to cut-out wind speeds.Blade elasticity is also included in the simulations.
The sampled time-series of blade strip AoA and relative velocity (,   )  are extracted from a set of time-dependent simulations.Each simulation considers the WT operating for 10 min in turbulent wind with a mean wind speed   at hub height ranging from 4 to 24 m/s with step of 1 m/s.Thus, 21 simulations are performed, with each simulation using a time-step of 0.01 s and storing the solution every 0.1 s.Therefore,   = 6000 samples are obtained for each mean wind speed.The turbulent wind for the inflow is generated using Turb-Sim [45], producing stochastic wind with mean vertical shear defined by Eq. ( 28) and speed fluctuations based on the Kaimal spectrum, using the turbulence intensity defined by IEC 61400 for Class I-B WTs.
The calculations of the LEE damage are conducted on the blade sections from 70% to 95% of the blade tip radius, as in [32].In this range, the blade is discretized with six strips, highlighted with different colors in Fig. 14(a), all having the cross-section defined by the NACA64 3 -618 airfoil, reported in Fig. 14(c).This figure also shows the airfoil geometry discretization adopted for the LEE damage calculation, which consists of 400 segments along the entire airfoil.Fig. 14(b) indicates the origin of the curvilinear coordinate ∕.

Computation of impact pattern
As highlighted in Eq. ( 1), the LEE damage depends on   and   .The value of   depends on the material properties and erosion model (discussed in Section 3.1), and the specific impact characteristics, i.e. droplet size, impact velocity and impact angle.The value of   is determined by the rain and wind conditions and the relative wind speed and AoA of the blade strip being considered.Therefore, the prediction of the damage incubation , the incubation time 1∕ and the growth of the erosion depth  after completion of the damage incubation, requires considering all possible impact conditions on the blade during the period of operation.
The joint FDFs of droplet size and wind speed were defined in Section 3.2.For given mean wind speed, however, the blade sections experience aerodynamic variations due to turbulence and only partly attenuated by the WT control.Thus, a method is required to use the rain and wind FDF   (  , ) for obtaining   , the number of droplets actually impacting the blade surface.
A simple method for estimating the minimum incubation period neglects the airfoil geometry, does not resolve blade aerodynamics and assumes a normal impact at the LE.In order to calculate the number of impacting droplets at a given blade radius , the FDF   (  , ) is multiplied by the ratio of local blade velocity and droplet terminal velocity.Applications of this type of method are reported in [5,31].The approach above can be refined by considering an impingement efficiency, namely a corrective factor accounting for the modification of the trajectory and impact characteristics of the rain droplets caused by aerodynamic forces [10,12].The impingement efficiency is calculated with the method proposed by Langmuir and Blodgett [22,23].This factor depends on the droplet size; in the complete erosion analysis, the impingement efficiency is evaluated for each droplet diameter, and used as a multiplicative factor to correct the impact velocity computed with the local blade speed and the droplet terminal velocity.
Approaches that capture more physical features of the interactions of hydrometeors and blade surface can also be used to determine the impact parameters.In this study, two different approaches, both considering the airfoil geometry of the blade strip, are adopted to compute the set of impact variables   ,  , and  for each segment  of the strip airfoil.Consideration of the strip geometry enables both accounting for the local airfoil curvature on the direction of the droplet impact and considering the dependence of coating durability and mass removal on the curvilinear position along the strip airfoil.The two approaches differ only because one approach neglects the impact of the aerodynamic forces on the rain droplets in the calculation of their trajectories and impact point, whereas the other approach considers also the aforementioned aerodynamic forces.
As indicated in the block diagram of Fig. 1, the characteristic parameters of the droplet impact on strip segment  for mean wind speed   are calculated for each AoA   and relative velocity  , at instant  of the time-series generated with the turbulent aerodynamic analysis (as discussed in Section 3.3).
Using the method which does not consider the impact of blade aerodynamics on the droplet trajectories, the impact angle is computed as   (∕) = arccos ( ) at each instant of the time-series, where

𝐕 𝐛
|  | = (cos , sin ) is the unit vector of the direction of the relative velocity, and (∕) is the unit normal vector of each strip segment .The impact velocity is taken to be  , .The variable    is obtained with the equation: Here, is a scaling factor accounting for the number of turbulent velocity samples for the strip aerodynamics   at mean wind   , and to convert from the droplet terminal speed   to the blade relative velocity by scaling the FDF   (  , ), as described also in [32].The modeling features above refer to the approach that does not include blade aerodynamics in the erosion assessment.A higher-fidelity improvement of this method also resolves the effects of the aerodynamic forces on the trajectories and impact characteristics of the Fig. 13.FDFs of wind speed and rain droplet diameter for assessing sensitivity of predicted coating lifetime to DSD source and FDF assembly.Same wind data from WT SCADA are used in all three cases.Map of Fig. 13(a) is equal to that of Fig. 11(a).
rain droplets [25].The first stage of this method requires computing a database of CFD and particle tracking simulations for each airfoil considered in the erosion analysis.The dataset of each airfoil consists of CFD analyses at different far field velocity and AoA, and, for each airfoil, a set of Lagrangian particle tracking simulations.All Lagrangian analyses use the same particle flow rate but different droplet diameter.The output of interest of each Lagrangian simulation is the distributions along the airfoil of the impact angle   , the impact velocity, the number of impacts per unit surface normalized by the number of droplets per unit surface injected in the simulation.For each airfoil, a surrogate data-driven model consisting of three transfer functions is then trained by using the database defined above.The surrogate model returns the same type of output of the underlying Lagrangian simulations for a flow of particles of diameter  for user-given relative velocity and AoA.The detailed methodology and verification of the metamodels are reported in [25].

Results
The LEE analyses and comparisons reported herein refer to realistic onshore WT operation, whereby the NREL 5 MW turbine discussed in Section 3.3 is assumed to be operational at the Hazelrigg weather station site.Some of the analyses below use part of or all the yearly wind speed and rain data available at the site, including a joint FDF of wind speed and rain droplet diameter based om wind speed data from the nacelle anemometer of a 2.3 MW operating at the site, and the DSD recorded by a laser disdrometer.
The section is made up of four subsections.The reference or baseline erosion analysis, which includes the assessment of both the coating durability (i.e.end of incubation time) and the time history of the mass removal of the coating, is presented in Section 4.1.This reference analysis illustrates the capabilities of the methods presented in this article, presenting and defining all variables and parameters used in the subsequent comparative analyses.The variability of the coating life prediction associated with different approaches to the calculation of the erosion strength and material weathering is analyzed in Section 4.2.The scatter of LEE predictions due to considering or neglecting the effect of the aerodynamic forces on the trajectories of the rain droplets impinging on the blade, and the source and combination of wind and rain data is analyzed in Section 4.3.Finally, the variability of the time to complete failure of the LEP arising from the choice of the modeling options and data types mentioned above is quantified in Section 4.4.

Reference erosion assessment
The reference erosion analysis uses the joint FDF of wind speeds and rain droplet sizes based on the nacelle anemometer of the Hazeltigg 2.3 MW WT and the DiVeN DSD.The aerodynamic effect on the trajectories of the rain droplets is included in the analysis by using the approach summarized in Section 3.4.The erosion characteristics of the coating are determined using the DNV-GL-RP-0573 guidelines reported in Section 3.1.2),an approach that was labeled Case 2.
Table 5 reports the predicted coating lifetime, namely the time at which the first erosion damage occurs, for the blade sections of the six outermost blade strips defined in Section 3.3.Variables   and  denote, respectively, mean distance of the strip from the rotor center and chord of the strip.The results highlight that the coating life is only around one year in the tip region, at 56 m <   < 59 m, where the blade speed is highest.For all six strips, the point at which the first damage occurs is on the pressure side, at distance ∕ = 0.0022 from the LE, corresponding to the first pressure side segment.The independence of this position on the strip, i.e. the radial position along the blade, occurs because the variability range of the mean AoA, similar for all six strips, is fairly limited below rated wind speed, where the majority of the rainfall occurred.
The maps of Fig. 15 report the fraction of the damage incubation (∕) after one year of operation computed using the measurementsbased joint FDF of wind speeds and rain droplet sizes for the six outermost strips.Each component (∕)  of the map gives the damage fraction accumulated at mean wind speed (  )  for droplets of diameter   .For each strip, the durability values reported in Table 5 are the inverse of the summation of (∕)  reported in the maps of Fig. 15 multiplied by the number of hours in a year.The pattern of the damage distribution in Fig. 15 is qualitatively similar for all strips, as a consequence of similar characteristics in terms of mean impact pattern over the surface.Indeed, after a given time of operation, the damage extent depends mainly on the relative velocity magnitude, which, in turn, depends on the radial position along the blade.
Fig. 16 presents the map (∕)∕  of strip 5, where   denotes the joint FDF of wind speeds and rain droplet sizes used as input of this erosion analysis.Each element [(∕)∕  ]  of the map is obtained by performing an element-by-element division of the (∕) and   FDFs, i.e. [(∕)∕  ]  = (∕)  ∕(  )  .The map displayed is a measure of the erosion potential depending on the droplet size, wind speed, material strength and WT rotor geometry and control.The largest values of (∕)∕  occur for the largest droplet diameters when the WT rotor speed achieves its maximum at   around 11.4 m/s.
Once the end of incubation is reached, erosion begins to remove material from the coating.This process is modeled using the approach outlined in Section 2. Since RET data were available only for mass removal from the coating, the computations below focus on the removal of the coating material, neglecting erosion of the substrate material.Therefore, the process of computing the mass removal ends once the erosion damage depth equals the thickness of the coating ℎ  = 600 μm.
Fig. 17 presents the profiles of erosion depth  along the airfoil of all six strips after one year of operation, with positive and negative values of ∕ referring to the pressure and suction side, respectively.Strips 1 to 4 are not affected by significant erosion, whereas strips 5 and 6 undergo erosion, with depth level differing significantly due to different blade speed.It is noted that, when considering the dimensional curvilinear coordinate , the actual curvilinear extent of LE erosion for strips 5 may become larger than for strip 6, because the ratio between the chords of strip 6 and 5 is about 2/3.
With regard to the implications of the computed erosion profiles on aerodynamic performance, it is noted that the maximum erosion depth per unit chord of strip 5 is 0.01%, whereas that of strip 6 is 0.04%, four times larger.Since the damage depth per unit chord and also the curvilinear extent of erosion per unit chord on the suction side strongly affect the aerodynamic performance of the blade and, thus, WT power [2,3], the analyses of the erosion characteristics of strips 5 and 6 highlight the importance of adequately protecting this blade from LEE.
Fig. 18 reports the computed damage progression of strip 5 over 30 years of operation.This analysis is accomplished by scaling the FDF   by a duration of 30 years, as done in previous analyses of this type [16].It is noted that after a certain time (about 10 years in this example), the damage does not grow significantly in the chordwise direction, particularly on the suction side, due to the probability of droplet normal impact decreasing moving away from the LE.

Sensitivity of incubation time to method for determining erosion strength, and material weathering
The analysis of this section aims at determining the impact of the method used to determine the erosion strength and material weathering on the prediction of the coating durability.In order to accomplish this, the four set-ups or Cases described in Section 3.1 are considered.The wind and rain data and the approach used to calculate the impact characteristics of the rain droplet are those of the reference erosion assessment for all set-ups.
Table 6 reports the LE coating durability of the six strips predicted using the four set-ups.Significant variations of the four predictions  are observed.Cross-comparing Case 1, which uses the fundamental material properties to obtain the coating erosion strength (Springer model), and Case 2, which used RET data according to the DNV-GL-RP-0573 guidelines to obtain the coating erosion strength, reveals that the former method significantly underpredicts the coating life.This result is consistent with the V-N curves of the two Cases reported in Fig. 7.The LE coating life predicted in Case 2 is between 1.8 to 2.2 times longer than in Case 1, with the mismatch between the two predictions growing with the blade radius.The cross-comparison of the data of Table 6 for Case 3, with weathered coating, and Case 4, without weathering, illustrates the predicted impact of weathering on the coating material lifetime.The coating not exposed to weathering has erosion life about 1.4 times longer than the weathered material.Equivalently, this result implies that in the considered experiment, weathering reduces the coating lifetime by nearly 30%.The reduction appears to be fairly independent of the radius, i.e. the blade peripheral speed.

Sensitivity of incubation time to blade aerodynamics, and wind and rain modeling
Using the set-up of the reference erosion assessment presented in Section 4.1 to determine the erosion characteristics of the coating, the comparative analyses presented below aim at assessing the impact of blade aerodynamics, wind data variability, and wind and rain data type and FDF generation on the predicted durability of the LE coating.
The considered analysis set-ups and the computed coating durability of each strip are summarized in Table 7.The first column indicates if aerodynamic effects are included in the calculation of the trajectories of the water droplets impinging on the blade.The second column indicates if the annual rain droplet DSD is based on the DiVeN measurements or the notional Best DSD.The third column indicates if the wind data used to generate the wind and rain FDF are based on the WT nacelle Comparing the results of the first two rows of Table 7 shows that the prediction of the coating life is slightly more conservative when the alterations of the droplet impact pattern caused by the aerodynamic deflection of the droplet trajectories are neglected.Neglecting the aerodynamic forces, the predicted coating life is about 12% shorter than when aerodynamic effects are considered for all blade strips, i.e. independently of the radial location.The shorter life predictions obtained without blade aerodynamics are due to the fact that the smaller droplets are deflected in the direction of the air streamlines, an effect which results in an increase of the impact angle and, thus, the normal component of the impact velocity.
Comparing the data of the first and third rows of Table 7 shows a much smaller variation of the predicted coating life using either an annual wind speed record based on the SCADA wind data or met mast wind data extrapolated to the hub height for generating the joint FDF   (  , ).The use of the turbine nacelle wind data rather than the wind data measured by the met mast and extrapolated at the nacelle height results in a reduction of about 3% of the coating life, fairly constant along the blade radius.It is remarkable that the significant difference between these two wind speed records observed in Fig. 10 yield a small difference in coating life prediction.
The results of the fourth row of Table 7 refer to the case in which the rain load and wind FDF is generated following the guidelines of the DNV-GL-RP-0573 standard.In this case, the only site-specific data for generating the FDF are the annual mean wind speed and shape factor of the Weibull distribution fitting the measured wind data, and the measured annual rainfall.The cross-comparison of this set-up and the reference erosion assessment highlights that the predicted coating life is 70% longer at all radii when using uncorrelated and partly site-agnostic wind and rain data.This conclusion, however, is sitedependent.Sites characterized by annual rainfalls similar to that at Hazelrigg, but featuring a significantly different DSD, may reverse this trend.The conclusion arrived at for the Hazelrigg site highlights the importance of using measured synchronous and co-located wind and rain measurements, whenever there are available, for the erosion analysis.The outcome of the comparative analysis above may be explained by considering the measured and modeled FDF of wind and rain data reported in Figs.13(a) and 13(b), respectively.The measured FDF shows that the peak area of the distribution occurs at winds between 3 and 12 m/s, a range in which the WT operated many hours per year, whereas the modeled FDF has the peak area shifted towards lower wind speeds, where the WT is shut down, or operates with low blade speeds.Furthermore, the measured FDF highlights a greater number of larger (more damaging) droplets than the modeled FDF, consistently with the DSDs of Fig. 12.
The analysis summary of the fifth row of Table 7 refers to a rain and wind FDF generated in accordance with the guidelines of the DNV-GL-RP-0573 standard, except for the fact that the Best DSD is replaced with the measured DiVen DSD.The predicted coating lives, nearly 90% longer than those predicted by the reference analysis, are comparable to those of the fourth row, which refers to the standard DNV analysis.This outcome indicates that, for the considered site, the large difference of coating life prediction using either a measured joint wind and rain FDF or the DNV standard, is due primarily to the latter approach attributing the same DSD to all considered wind speeds and neglecting the timecorrelation of wind and rain data.Conversely, the differences between the DiVeN and Best DSDs visible in Fig. 12 appear to play a minor role in the coating life prediction.

Comparative assessment of LEP material loss
Tables 6 and 7 compared the incubation times of the six outboard strips obtained by using all considered material erosion models, coating states, and wind and rain data source.This section assesses the sensitivity of the erosion damage, i.e. the extent of the removed coating mass, to the use of the same LEE analysis set-ups.The summary of these comparative analysis is presented in Table 8, which shows the extent of coating removal at the LE predicted by all analysis set-ups after two years of operation.Each analysis is performed by simply doubling the bin content of the wind and rain FDF corresponding to the particular set-up.The erosion extent along the LE is defined by the positions at which material removal starts and ends using the normalized curvilinear coordinates ∕, with positive ∕ denoting a point on the pressure side and negative ∕ denoting a point on the suction side.Starting from the trailing edge and moving around one side of the airfoil to reach the trailing edge from the opposite side, the beginning of the erosion patch corresponds to the first point with nonzero  and the end corresponds to the last point with nonzero .
The results of Table 8 feature the same trends observed for the variability of the incubation time with the analysis set-up highlighted in Tables 6 and 7.For instance, the erosion model affects significantly also the coating's mass loss, since the use of Springer model (first row of the table) results in larger mass removal than the use of the DNV guidelines-based material erosion model (second line of the table), in line with the comparison of the incubation time in the two cases.Similarly, the mass loss recorded by using synchronous collocated measured data of wind speed and rain droplet size (second column of the table) is notably larger than determined by using the notional model of the industrial recommended practice (seventh row of the table).This is the same trend observed for the incubation time in the two cases, as seen in Table 7.
The comparison of the results of the analysis set-ups with and without aerodynamic forces (second and fifth rows of Table 6, respectively), show an interesting trend: for the outermost strips (strips 5 and 6) the chordwise extent on erosion is larger when blade aerodynamics is included in the analysis.This is due to the deflection of the trajectories of the smaller rain droplets approaching the LE in the region above the stagnation streamline.This effect results in these droplet hitting the suction side further downstream of the LE with respect to when blade aerodynamics is not included.This phenomenon does not occur on the pressure side because local airflow accelerations are notably smaller.The phenomenon described above does not occur on the suction side of the other strips because the momentum of the smaller droplets is not sufficiently large to promote erosion.Despite the fact that the overall curvilinear extent of erosion is larger when aerodynamic forces are included, the mass loss is lower.This is because erosion is shallower when aerodynamics is included due to the aerodynamic deceleration of the droplets and the deviation of their impact angle from 90 • .The resulting smaller erosion depth outweighs the larger curvilinear extent of erosion, and the mass loss is therefore smaller than in the analysis that does not include blade aerodynamics.

Turbine AEP loss due to observed LEE
In order to estimate the AEP losses associated with the LEE levels predicted by the reference erosion analyses, the AEP loss prediction system (ALPS) framework [1] has been used.Three time-consecutive damage configurations have been analyzed.The first stage of the geometric damage considers the LE coating material removed after one year of operation.The curvilinear abscissas defining the boundaries of the removed material, not reported herein for brevity, are calculated after one year of operation.Since the erosion analyses cap the maximum erosion depth to only 0.6 mm, the equivalent sand grain roughness model [46] is used in ALPS to account for the aerodynamic performance reduction of the eroded blade sections.
The second erosion stage considers the effective LE geometry after two years of operation.The boundaries of the removed mass areas are defined by the curvilinear abscissas in the second raw of Table 8.The erosion depth is still assumed to be 0.6 mm, and the equivalent sand grain roughness model is used also in this case.For both damages above, the actual roughness  is taken to be 0.6 mm, whereas the equivalent sand grain roughness   corresponds to a ratio   ∕ = 5.This choice of this ratio is made both because it is close to the mean value of the range considered in many studies and applications [47], and because it is a level at which the WT performance reduction does no longer depend significantly on   ∕ [2].
Since the assumption that the erosion depth remains at 0.6 mm for relatively long times after the LE coating has been completely removed, a third more realistic damage scenario is considered for the second year, namely one defined by the same curvilinear coordinates of the second damage, but an erosion depth of 3 mm.This corresponds to considering initial erosion of the composite substrate, which is typically much less resistant to erosion than the coating.The third damage level is analyzed by resolving its geometry in the aerodynamic analysis making use of a chordwise LE groove representation [1].The curves of lift and drag coefficients of the nominal and eroded variants of the NACA64 3 -618 airfoil for the three ALPS analyses above are determined using the database reported in [48].
The AEP estimates of the considered 5 MW WT also require defining the turbulence intensity (TI) for each mean wind speed.Similarly to the AEP analyses of [2], the TI data are based on the International Standard IEC61400 standard.Assuming an onshore wind of class B, the TI for each wind speed is given by   =    (0.75 + 5.6∕  ) where    = 14%.The AEP losses associated with the three considered damage stages are referred to the AEP of 10.852 GWh of the nominal WT, and are reported in Table 9.These data show that at the end of year 1 a modest AEP loss of only about 0.1% is recorded, due to the fact that material removal has affected by a very small amount only the two outermost blade strips.At the end of year 2, conversely, the AEP loss estimate varies between 0.4 and 1.4%, with the variability being due on the uncertainty depth of the LE substrate.

Conclusions
The article presented novel research on the prediction technology of rain erosion of WT blades.The developed multidisciplinary and multiphysics LEE analysis framework can use fundamental material properties or RET data for assessing LEE of real WTs, and allows considering the effects of site-specific wind and rain climates, blade geometry and blade aerodynamics on the erosion assessment.The impact of several modeling and data choices on the predictions was investigated.The geographical location of all analyses was a Northwest England onshore site for which comprehensive wind and rain data are available.All analyses were performed for a 5 MW WT, since the WT geometry and control data of a 2.3 MW WT operating at the considered site were not available.
Sensitivity analyses of the impact of modeling choices on LEE predictions considered the material erosion model and the use of blade aerodynamics for better resolving the impact characteristics of rain droplets.Parametric analyses of the impact of data type choices considered the definition of wind and rain data, based either on comprehensive site measurements or data models generated in accordance with sector standards, the source of wind speed data, using recordings of the WT nacelle anemometer or readings extrapolated from a shorted met mast, and the impact of coating material weathering.
The reference LEE analysis used annual wind speed data measured at 90 m above the ground by the nacelle anemometer of a 2.3 MW WT, and synchronous disdrometer data at the same site, used to generate a joint FDF of wind speed and rain droplet size.The reference analysis considered a typical industrial coating material of known material properties, used the RET-based erosion model of the DNV standard, and included the effect of blade aerodynamics in the prediction of the droplet trajectory and impact patterns.
The reference LEE analysis showed that the incubation time of the first LE damage, at the outboard part of the blade, occurred after 5509 h of operations, i.e. about 7 month and 3 weeks.The alternative analysis using the fundamental material properties predicted an incubation time of about half that of the reference analysis at the same blade radius.This underestimation of the incubation time is consistent with the findings of previous RET data-and Springer model-based comparative studies.It is noted, however, that applying the two erosion prediction methods to a WT analysis including real operation effects (e.g.droplet size variability and blade aerodynamics), the aforementioned ratio of about 2 between predicted incubation times of the WT blade decreased with respect to that of about 3 determined for the controlled environment of the swirling arm RET rig.
Modifying the reference analysis by neglecting the effects of blade aerodynamics on the droplet impact pattern resulted in an erosion life reduction of about 12%.
The sensitivity of the predicted erosion life to the wind and rain model and data was found to be significant, since the erosion analysis performed by using the wind and rain model recommended by the sector standard predicted an erosion life equal to 1.7 times that predicted by the reference analysis using the measurement-based joint FDF of wind and rain data of the considered site.This disparity was found to be due mostly to the lack of time-correlation between wind speed and droplet size data in the recommended practice, rather than the differences between the measured and notional DSDs.The use of a joint wind and rain FDF based on synchronous co-located anemometer and disdrometer measurements also enabled the identification of weather conditions critical for LEE.
The life prediction estimate using the nacelle or extrapolated 10-m wind speed data differed by only 2%.Material weathering due to UV exposure during WT operation was found to reduce the LEE coating life by about 30% with respect to the reference analysis, which did not include weathering.
The AEP loss estimate corresponding to the LEE damages predicted by the reference erosion analysis after two years of operation varies between 0.4 and 1.4%.The relatively wide variability range ensues from the uncertainty on erosion depth, capped at the coating thickness in the performed erosion analyses.It is envisaged, however, that erosion of the substrate, expected to progress at a higher rate than that of the coating, results in the AEP loss being closer to the upper end of the indicated range.
The presented study identifies several large sources of uncertainty in a state-of-the-art rain LEE prediction framework, thus providing guidelines on improvements in terms of required fidelity and input data type of each module.Discussed results highlight several future research and development avenues.Knowledge of the critical erosion conditions enabled by using site-specific and measurement-based joint FDF of wind and rain data is paramount for developing further and, ultimately, deploying erosion mitigation strategies, such as the socalled 'erosion-safe WT operation mode' proposed in [29].Presented results also highlight the risk of mispredicting LE coating life and material removal rates by using notional rather than measurementbased joint FDF of wind and rain data.Although the misprediction reported in this study is due to the lack of time-correlation of rain and wind data in the notional model of the recommended practice, it is envisaged that, for other sites, significant prediction errors may occur due to large differences between the notional DSD and the sitespecific DSD.This emphasizes the need for developing atlases of erosion potential, and deploying the experimental and simulation technology for characterizing atmospheric precipitations.

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.

A
.Castorrini et al.

Fig. 2 .
Fig. 2. Multilayer LEP system structure obtained with post-mold application, and RET sample configuration based on DNV-GL-RP-0171 guidelines.

Fig. 3 .
Fig. 3.One specimen without weathering at initial (left) and final (right) time of test set 1 A ( = 1159 RPM).Reported length corresponds to speed range 130-140 m/s.

Fig. 9 .
Fig. 9. Rainfall rate frequency distributions generated by using the TB data and by integrating the DiVeN droplet data over 1-h intervals.

Fig. 10 .
Fig. 10.Weibull PDF based on WT SCADA data and extrapolated met mast measurements.

Fig. 11 .Fig. 12 .
Fig. 11.Joint FDFs of wind speed and rain droplet diameter using DiVeN rain measurements with wind speed data based on WT SCADA data and extrapolated met mast measurements.

Fig. 15 .
Fig. 15.Maps of damage incubation fraction for six outermost blade strips after one year of operation predicted by reference erosion analysis.

Fig. 18 .
Fig. 18.Erosion depth progression of strip 5 over 30 years of operation predicted by reference erosion analysis.

Table 1
Reference RET operational parameters.

Table 3
LEP and water droplet material properties.

Table 4
Values of erosion strength and exponent  for Cases 1-4 with  = 2.6 mm in RET experiment.[GPa][GPa]( Fig. 7. V-N RET test sets and fitted curves for Cases 1-4.

Table 5
Coating lifetime prediction of considered strips for reference erosion analysis.

Table 6
Coating lifetime prediction of considered strips for different erosion and material models.

Table 7
Coating lifetime prediction of considered strips for different approaches to blade aerodynamics and wind and rain FDF generation.

Table 8
Prediction of LEE curvilinear extent of considered strips after two years for all erosion analysis set-ups.SCADA) or the met mast records extrapolated at hub height.The information on the construction of the wind and rain FDF is also contained in the second column.The isolated label 'DiVeN' indicates the use of a joint FDF generated from co-located synchronous wind speed and droplet size measurements, whereas the label 'DNV' indicates the use of an FDF generated in accordance with the DNV-GL-RP-0573 guidelines.Label 'DNV(Best)' refers to the case in which the site-agnostic Best DSD is used, which corresponds to entirely adhering to the DNV recommendations.Label 'DNV (DiVen)' refers, instead, to the case in which the Best DSD is replaced by the measured DiVen DSD, but the same DSD is attributed to all considered wind speeds, as recommended by the DNV standard.

Table 9
Estimates of AEP losses in first two years for LEE damages based on reference erosion analysis.Labels 'mod.' and 'grv.' refer to modeled roughness and geometrically resolved LE grooves, respectively.Modeled roughness analyses use   ∕ = 5.