We propose a data-driven model to predict the short-term fatigue Damage Equivalent Loads (DEL) on a wake-a ected wind turbine based on wind eld in ow sensors and/or loads sensors deployed on an adjacent up-wind wind turbine. Gaussian Process Regression (GPR) with Bayesian hyperpa- rameters calibration is proposed to obtain a surrogate from input random variables to output DELs in the blades and towers of the up-wind and wake-a ected wind turbines. A sensitivity analysis based on the hyperparameters of the GPR and Kullback-Leibler divergence is conducted to assess the e ect of di erent input on the obtained DELs. We provide qualitative recommendations for a minimal set of necessary and sucient input random variables to minimize the error in the DEL predictions on the wake-a ected wind turbine. Extensive simulations are performed comprising di erent random variables, including wind speed, turbulence intensity, shear exponent and in ow horizontal skewness. Furthermore, we include random variables related to the blades lift and drag coecients with direct impact on the rotor aerodynamic induction, which governs the evolution and transport of the meandering wake. In addition, di erent spacing between the wind turbines and W ohler exponents for calculation of DELs are considered. The maximum prediction normalized mean squared error, obtained in the tower base DELs in the fore-aft direction of the wake a ected wind turbine, is less than 4%. In the case of the blade root DELs, the overall prediction error is less than 1%. The proposed scheme promotes utilization of sparse structural monitoring (loads) measurements for improving diagnostics on wake-a ected turbines. Keywords: wind turbine, fatigue, wake, uncertainty, Gaussian Process Regression, Bayesian calibration, sensitivity analysis, virtual sensing Nomenclature Wind shear exponent E Expected value of a random variable V Variance of a random variable D Rotor diameter In ow horizontal skewness  Turbulence CD Aerodynamic drag coecient CL Aerodynamic lift coecient DEL Damage Equivalent Load GPR Gaussian Process Regression mBWh oler exponent for blade (composites) mTWh oler exponent for tower (welded steel) ML Machine Learning RV Random Variable 1Authors made equal contributions to the study and the publication Preprint submitted to engrXiv November 5, 2020 Ti Turbulence intensity U Mean wind speed WT Wind Turbine 2 1. Introduction The load e ects and degradation of structural components of wind turbines (WTs) do not uniformly evolve across a wind farm due to, largely, wake e ects and variability in the in ow conditions (and waves for o shore wind farms [25]). The assessment of fatigue damage accumulation on assumption of availability of direct load e ects measurements on all main structural components across all wind turbines within a wind farm is not realistic. The assumption of availability of high- delity aero-servo-elastic simulators of the investigated wind turbines coupled to site-speci c in ow and wake models is convenient, but often not borne out of the actual reality experienced by wind farm owners and operators. At best, the designer/manufacturer of the wind turbine might make a single comprehensive one-o set of output load e ects of such wind farm simulations available to the owner/operator of the wind farm. A further and, from the view of this paper, perhaps more important hindrance lies in the precise estimation of structural response signals that are typically not available in the standard data emanating from the Supervisory Control and Data Acquisition (SCADA) monitoring systems embedded on wind turbines. A number of works therefore attempt condition monitoring on the basis of SCADA availability [47, 26]. However, when structural monitoring information becomes available [16], then this can be exploited as a more direct proxy to diagnose sudden damage [8], or to further accurately assess the remaining useful lifetime of a wind turbine in a given wind farm [17, 19, 33, 44]. We therefore posit that considerable improvements to the operation, maintenance and prediction of remaining useful life of a wind turbine can be accomplished by delivering access to transparent, simple, yet powerful and interpretable data-driven predictive models. Such models could be trained, tuned and updated via (fairly) easily accessible and cheap structural response observations from a limited number of appropriately selected wind turbines in a given wind farm [31]. The work presented herein focuses on the e ects induced by wakes on the Damage Equivalent Loads (DELs) on structural components of WTs. The fatigue load variation within an o shore wind farm is primarily a product of wake-induced ow disturbances [42]. The form and severity of the wake-induced turbulence and de cit depend on multiple factors, which are in this work grouped into three classes:(a) ambient conditions (wind speed, ambient turbulence, atmospheric stability, wave height, and others);(b) wind turbine operational regime (rotor thrust, rotational speed, and power set point); and (c) relative position of the wake source(s) with respect to the disturbed turbine. Recent literature has attempted to tackle the fatigue variability issue, by delivering predictive frameworks that capitalize on availability of SCADA data. The approaches in delivering such predictive models may be distinguished in terms of two main categories, namely the physics-based and data-driven classes. Initiating from a physics-based approach, research in [12] and in similar works [41, 37] exploits a surrogate approach, relying on Polynomial Chaos Expansions (PCE) and Arti cial Neural networks (ANN), trained on pre-simulated load scenarios, to predict the fatigue load variation on WTs for a wind farm with arbitrary layout under wake e ects. ANNs are shown to outperform PCE in terms of prediction accuracy and computational speed [40], albeit being prone to over tting, while further require signi cantly more data for achieving acceptable performance. Both methods allow for obtaining analytical derivatives, which is a useful trait in optimization and sensitivity analysis. In [11] the performance of ve surrogate models is assessed by comparing site-speci c lifetime fatigue load predictions at 10 sites using an aeroelastic model of an individual DTU 10MW reference wind turbine. The compared methods include PCE, quadratic response surface, universal Kriging, importance sampling, and nearest-neighbor interpolation. The authors argue that PCE-based (and Kriging) models may sometimes have a practical advantage over ANNs, due to the \white-box" features { such as being able to track separate contributions to variance (and uncertainty). Research in [15] proposed a procedure for producing a lifetime fatigue load variation map within an o shore wind farm. Factors such as ten-minute average free wind speed, free wind direction, ambient turbulence, farm layout and wake e ects, wave height, peak period, and alignment with the wind were considered. The procedure relies on direct aero-elastic simulations of the whole wind farm including wake e ects using the DWM model. A similar approach was adopted in the work of Tagliatti [43]. This mapping is not directly extendable for use with continuous and long-term structural health monitoring data (SHM). In a purely data-driven scheme, which does not take aeroelastic analyses into account, Pap- atheou et al. [34] focus on power prediction for the Lillgrund wind farm [2] for the purpose of condition monitoring and fault detection. They adopt both ANNs and Gaussian Processes (GPs) for producing individual and population-based power curves. They then attempt prediction of the power produced on individual WTs based on measurements extracted from other turbines in the farm. A comparison between neural networks and GPs reveals no signi cant di erence in terms of precision, but showcases the inherent ability of the GPs to produce probabilistic bounds. Woo 3 et al. [46] propose a Multi-Tasks Convolutional Long Short-Term Memory Network approach to simultaneously predict the energy output and structural load from the target wind turbine, while modeling the spatio-temporal structure of the input wind ow. The work is veri ed on simulations from a stand-alone NREL-5MW onshore reference wind turbine. The predictions are delivered in a short-term horizon, i.e., few seconds ahead. In both [34, 46] wake e ects are not considered. Wake is tackled in [32], where a trained Variational Autoencoder (VAE) is exploited to map the high dimensional correlated stochastic variables over the wind-farm, such as power production and wind speed, to a parametric probability distribution of much lower dimensionality, with the ultimate goal of condition monitoring. In a method that attempts to fuse physics-based simulators with data, in what concerns the train- ing of predictive ML models, Dimitrov et al. [13] propose to use a combination of limited SCADA based measurements and wind turbine/farm simulations. An arti cial intelligence framework is trained to forecast the future performance of the wind turbine and the fatigue life consumption of its components. If SCADA measurements, such as measured power production, wind speed and rotor speed, are available, the load mapping can be realized by training a data-driven regression model using e.g. arti cial neural networks (ANN). In the absence of actual loads data, an aeroelastic model of the turbine can be used to generate a synthetic data set to serve for training. A com- parison of the normalized damage equivalent blade root ap moments between measurements and simulations shows a discrepancy of the order of 10% 15% for some of the operating wind speeds. Park & Park [35] present an attempt to fuse data with engineering principles via a physics-induced graph neural network (PGNN) model able to estimate the power outputs of all wind turbines in any layout under any wind conditions. An engineering wake interaction model serves as a basis function, which e ectively imposes physics-induced bias for modelling the interaction among wind turbines into the network structure. To clearly understand the role of the physics-induced weight function, the authors compare the PGNN performances to a purely data-induced approach, termed data-induced GNN (DGNN). When a target turbine in a wind farm experiences more complex wake interaction, the DGNN tends to overestimate power generations. A drawback of the proposed method is the inability to produce probabilistic output. The overview of existing literature reveals that, on the one hand, deep learning (DL) neural network based methods (Recurrent, Convolution, GraphNets, Auto-encoders, etc.) are suited to the problem at hand, but require special understanding and tuning of the network layers to reach adequate predictive results. Moreover, these often require large training datasets in order to avoid over- tting. Despite signi cant advancements in this eld, DL methods remain black-box in their very nature and comprehensive interpretable results remain dicult (for now). A number of surro- gate approaches are proposed but are focused on o ering an emulator of the load responses across wind turbines in a wind farm based on pre-compiled aeroelastic simulations. A lack is noted with respect to comprehensive e orts for predicting wake-induced load e ects on WTs. An important goal would be to predict the loads on unmeasured wake-a ected WTs via direct input loads and in ow measurements extracted from wake-free wind turbines in the same farm. This problem is particularly relevant in a realistic SHM setting for wind farms, since an ecient monitoring scheme, should involve carefully planned and, to the degree possible, sparse structural measurements. In this work, we propose adoption of classical Gaussian process Regression (GPR) with Bayesian learning of hyper-parameters. We argue that GPR o ers a number of advantages over deep learning, or some of the further surrogate modelling approaches, in that it delivers an elegant mathematical formulation and exact inference, it o ers a exible encoding of linear constraints, it has proven robust in small low dimensional input spaces and scalar univariate (non-time series) output [6, 7]. Perhaps the major advantage lies in the built-in feature for uncertainty quanti cation that enables e ective policies for data acquisition and experimental designs, including Bayesian approaches for hyperparameters optimization. On the downside, some limitations, which ought to be acknowledged include limited scalability to large data-sets and high dimensions, as well as limited expressivity and robustness to prior assumptions, especially in relation to the choice of kernels. The former consideration does not pose an issue for the analysis presented herein, as we do not deal with output time series datasets, but instead treat aggregated features (such as DELs). Furthermore, our input dimensional space is generally limited to few essential in ow and turbine response Random Variables (RVs). Regarding the second consideration on expressivity, the choice of GPR kernels is indeed a source of uncertainty. This can be tackled when kernels are introduced as a RV, as done in [3], or alternatively Kernel selection could be performed using Approximate Bayesian Computation, as done in [5]. In order to verify the proposed GP-based approach, an exploratory and therefore simpli ed analysis is adopted in this work, featuring an essential setup comprised of two interacting wind turbines; a rst WT situated up-wind, with the second positioned directly in the wake of the rst. This simple setup allows us to illustrate our ndings on utilization of the proposed framework by 4 means of easily interpretable results. Furthermore, we hypothesise that this setup is suitable for a signi cant number of small to medium onshore wind farms, such as for instance the layout shown in Figure 1. The wind rose indicates a narrow band of wind direction between North-North-East and North-North-West, resulting in single meandering wake eld amongst the turbines, which is the setup adopted in this paper. Figure 1: Layout of an onshore wind farm in complex terrain, located in central Greece. The main contributions of this work pertain to i) identi cation of dynamic di erences between up-wind and wake-a ected turbines; ii) development of GPR-based framework to estimate short- term fatigue DELs for the wake-a ected wind turbine, based o loads and in ow measurements on the up-wind turbine, and iii) recommendations for a minimal set of necessary and sucient input random variables to predict the short-term fatigue damage equivalent loads (DEL) on a wake- a ected wind turbine. In appropriately accounting for inherent uncertainties, beyond in ow and fatigue related uncertainties, in our design of experiments we inject direct aerodynamic uncertainties on the lift and drag coecients of the airfoil sections along the span of the blade, thus a ecting the rotor induction and consequently the dynamic wake evolution and transport. The remainder of this article is organized as follows. In Section 2 we describe the uncertainties and simulations setups. In section 3 we provide an interpretation of the main output from the numerical simulations especially with respect to the e ects of the uncertainties (Random Variables) in relation to the short-term DEL of various structural components on the wind turbines. In section 4 we elaborate on the framework of virtual fatigue diagnostics of the wake-a ected wind turbine via Gaussian Process Regression (GPR) model, and present the ensuing results in section 5. We nish with concluding discussions and outlook in section 6. 2. Uncertainty modeling and wake simulations setup In this section, we detail our uncertainty framework and the wake meandering aero-elastic simu- lations setup. Three categories of RVs are considered, namely: wind in ow RVs, aerodynamic RVs and fatigue RVs. 2.1. In ow RVs and their stochastic models The variation in the structural dynamic response of wind turbines is signi cantly dependent on the turbulent in ow wind eld conditions, including the mean wind speed, turbulence, wind shear, and in ow skewness. In accounting for these in uences, we introduce the following RVs in the simulations setup: Mean wind speed, U, turbulence intensity, Ti, wind shear, , and horizontal in ow skewness, . 5 The mean wind speed follows a Weibull distribution UWBL (AU;KU), truncated to [4  25]m=s, with parameters speci ed as follows: E(U) = 8:5; whereA U=2E(U)p KU= 2:0(1) The conditional dependence between the turbulence Uand the mean wind speed Uis de ned in the Normal Turbulence Model described in the wind turbine design standard [1]. Here, we elect to use a reference ambient turbulence intensity Iref=:16 (the expected value of the turbulence intensity at 15 m/s is called Iref). This dependency is given by the local statistical moments of ULN U;2 U as: E UjU =Iref(0:75u+ 3:8) V UjU = (1:4Iref)2(2) The wind pro le above ground level is expressed using the power law relationship, which de nes the mean wind speed Uat heightZabove ground as a function of the mean wind speed at hub heightUhmeasured at hub height Zhas reference: U Uh=Z Zh (3) where is a constant called the shear exponent. The conditional dependence between the wind shear exponent N ( ;2 ) and the mean wind speed Uis given by [14]: E jU = 0:088 (ln(u)1) V jU =1 u2 (4) We de ne a custom conditional dependence between the in ow horizontal skewness and the mean wind speed U, truncated to [11;11]deg:, N ( ;2 ): E jU =ln(u)3 V jU =15 u2 (5) 2.2. Aerodynamic RVs and their stochastic models A well know result from Blade Element Moment theory (BEM) links the aerodnynamic lift ( CL) and drag (CD) coecients to axial ( a) and tangential ( a0) induction factors as follows: a 1a=0(CLcos +CDsin ) 4 sin2 a0 1a=0(CLsin CDcos ) 4rsin2 (6) where0is the rotor solidity,  is the angle of the incoming relative wind with the rotor plane, tip speed ratio =!r V0,!is the rotor speed, ris the radial distance from the rotor center, and V0 is the freestream wind speed. The distribution of aerodynamic axial and tangential induction over the rotor essentially governs the evolution and transport of the wake [24]. Hence, our approach to a ecting aerodynamic induction is by introducing a stochastic model of the lift and drag co- ecients curves in BEM. Several sources of uncertainties a ect the lift and drag coecients with direct impact on aerodynamic induction. These uncertainties are associated with assessment of air- foil characteristics in wind tunnels, uncertainties due to 3D ow correction, uncertainties stemming from surface roughness, uncertainties related to the blade geometric distortions in manufacturing and handling, uncertainties related to the blade geometric distortions when de ected under load, uncertainties due to the e ects of Reynolds number, uncertainties associated with extending airfoil aerodynamic characteristics to post stall, and nally uncertainties stemming from the validation of airfoil data by eld full scale measurements. It is not possible to quantify the joint distribution of all these RVs and, as a result, a simpli ed approach is chosen via a stochastic model, as proposed in [4]. The stochastic model consists in parameterizing the lift coecient curve by the slope in the linear range@CL @ , the point indicating the start of the trailing edge separation ( AoATES;CL;TES ), the point of maximum lift ( AoAmax;CL;max ) and the point where the stall recovery is initiated 6 (AoASR;CL;SR). The drag coecient is several orders of magnitude smaller than the lift coecient for small angles of attack (below stall) and, thus, its impact is limited. Furthermore, it generally displays minor variability regardless of the airfoil type, geometry, or thickness to chord ratio. Conse- quently, the drag coecient is only parameterized by the point where minimum and maximum drag coecient occurs at AoA = 0andAoA =90, respectively. According to [4] the probabilistic distributions, expected values, coecient of variations and correlation coecients are assigned to the aforementioned parameters (for brevity, we do not repeat these here). Figure 2 shows samples of the stochastic lift and drag coecients curves. These perturbations result in modi ed CLand CDcurves that maintain the primary characteristics of the original aerodynamic polars, but dif- fer in both magnitude and feature location. Note that these synthetic aerodynamic lift and drag coecients curves are sampled independently from the wind in ow RVs. Furthermore, we vary the spacing between the up-wind and wake-a ected WTs, as shown in Table 1. Table 1: Spacing between wind turbines. Random variableDescription Probability DistributionParameters D Spacing in multiples of rotor diameters between turbinesDiscreteD= [3;5;8;11] -200 -100 0 100 200 Angle of Attack [deg.]-1.5-1-0.500.511.5Lift Coefficient [-]Generated Original (a)CL -15 -10 -5 0 5 10 15 Angle of Attack [deg.]0.0050.010.0150.020.0250.030.0350.040.0450.05Drag Coefficient [-]Generated Original (b)CD Figure 2: Samples of the stochastic lift and drag coecients CLandCDcurves of airfoil NACA 64-618 2.3. Sampling of the random variables We choose to sample the wind in ow and the aerodynamic RVs using the Sobol Quasi-Random sequences, which are designed to generate a sample that is uniformly distributed over the unit hypercube, i.e., as uniformly as possible over the multi-dimensional input space [39]. In total we sample 2048 Joint samples of the wind in ow RVs, as shown in Figure 3. For a given spacing between the up-wind and wake-a ected WTs, a sample of U,U, and , combined with a sample of stochastic CLandCD, we generate a realization of an in ow turbulent wind eld time series as input to the FAST-DWM aero-servo-elastic environment to simulate the corresponding aero-elastic response of the wind turbines structure. 2.4. Fatigue RVs and their stochastic models In this paper, we represent fatigue using the short-term fatigue damage equivalent load (DEL) concept. The advantage of the DEL is that it reduces a long history of random loads to one number, which makes it convenient to compare various load and operating scenarios [45]. In our probabilistic calculations the exponent of the S-N curve m(Wh oler exponent) is considered to be a RV. To get an impression of the in uence of the Wh oler exponent we compute the DEL based on a range of discrete Wh oler exponents for a given material as shown in Table 2. We assume that the blades composites Wh oler exponent varies between 9-13. We assume that the tower structural/welded steel Wh oler exponent varies between 3-4. The short-term fatigue damage equivalent loads (DEL) follow from the computed 10-min output time series response of the wind turbine: DEL = 1 NeqX ini(Di)m!1=m (7) 7 Figure 3: Samples from the joint wind in ow random variables. wherenidesignates the number of load cycles with range Diin theithrange bin of the fatigue load spectrum, and Neqspeci es the equivalent number of load cycles. Table 2: Fatigue related random variables. Random variableDescription Probability DistributionParameters mT Wh oler exponent for tower (welded steel)Discrete mT= [3;4] mB Wh oler exponent for blades (composites)Discrete mB= [9;10;11;12;13] 2.5. Dynamic wake meandering and aero-servo-elastic simulations setup Our dynamic wake meandering and aero-servo-elastic simulations setup is based on the coupled DWM [30] and FAST numerical models [20]. The simulations in DWM-FAST considered two reference NREL three-bladed up-wind, horizontal-axis WT [21] with 126 mrotor diameter, 5 MW rated power and hub height of 90 m. The rated power of 5 MW occurs at a wind speed of 11 :4m=s and a rotor speed of 12 :1RPM . We list some of the more important properties of the simulated wind turbine in Table 3. In the DWM setup, the two turbines are aligned as shown in Figure 4. Table 3: Properties of the NREL 5-MW reference wind turbine. Number of blades 3 Rotor diameter 126 m Hub height 90 m Rated power 5 MW Cut-in wind speed 3 m=s Cut-out wind speed 25 m=s Control Variable Speed, Collective Pitch Variable speed from cut-in to cut-out wind speed Variable pitch from cut-in to cut-out wind speed Rated wind speed 11 :4m=s Cut-in and rated RPM 6 :912:1RPM FAST is a wind-turbine-speci c time domain aeroelastic computer simulator that employs a combined modal and multibody dynamics formulation, adopting limited degrees of freedom (DOF). Since FAST models exible elements using a modal representation, the reliability of this represen- tation depends on the generation of accurate mode shapes by the engineer, which are then used 8 as input into FAST. Large structural elements, such as blades and tower models, are characterized by properties such as sti ness and mass per unit length to represent the exibility characteristics. FAST models the turbine using 24 DOF, including two blade- ap modes and one blade-edge mode per blade, two fore-aft and two side-to-side tower bending modes, nacelle yaw, the generator azimuth angle and the compliance in the drive train between the generator and hub/rotor. The aerodynamic model is based on the Unsteady Blade Element Momentum theory, including skew in ow, dynamic stall and generalized dynamic wake [10]. The Blade aerodynamic pro les' properties are provided as a-priori input and are used as lookup tables or for interpolation. The stochastic input wind eld uses the Kaimal turbulence model [22]. Aeroelastic simulations of WTs are stochastic largely due to the stochastic nature of the input wind eld; It is thus a common practice to generate a signi cant number of stochastic simulations for various operating and environmental conditions in order to cover variability on aeroelastic fatigue and extreme load analysis. Wind turbines located in wind farms experience a wind eld that is modi ed compared to the undisturbed ambient wind eld. A wake is characterized by a decrease in the mean wind speed and increase in wind speed uctuations (turbulence) behind a turbine. The downstream transport of a wake follows a stochastic pattern known as wake meandering (oscillations). It appears as an intermittent phenomenon, where winds at down-wind positions may be undisturbed for part of the time, but interrupted by episodes of intense turbulence and reduced mean speed as the wake hits the observation point [28]. Thus, a correct wind turbine load prediction requires the inclusion of the downstream evolution of wake de cit, the increased small-scale wake turbulence and the wake meandering. In this paper we choose to use the DWM wake model coupled to FAST following the NREL implementation. This coupling is well documented in [18]. The dynamic wake meander model coupled into FAST is used to model the up-wind (Turbine 1 in Figure 4) wind turbine's wake e ect on the structural dynamics of the down-wind wake-a ected turbine (Turbine 2 in Figure 4). While FAST is simulating an up-wind turbine, DWM calculates the wake de cit velocity, the meandered wake center positions with re- spect to time, and the added turbulence intensity due to the presence of the wake mixing. While a down-wind wake-a ected turbine is being simulated in FAST, the in ow wind to this wake-a ected turbine is modi ed based on the wake modelling results of its up-wind turbines. Thus, the e ect of the wakes can be re ected on the wake-a ected turbine according to its immediate wake [18]. It should be noted that the wake-induced load e ects are obtained under the assumptions underlying the DWM model, i.e, that the wake de cit behaves as a passive tracer following the transverse wind uctuations and that the ambient turbulence causing the meandering can be described by a Gaussian random turbulence model, such as the Mann model. Consequently, the FAST and DWM models might su er from model-form diciencies and lack of inclusion of some physics, which should not distract from the main objective and thrust of this work. Our simulations setup is limited to only two turbines, with ow down a row. We would expect to see more di erences in larger wind farms because of blockage and deep array e ects. u Ambient wind ﬁeld Aeroelas�c wake-free turbine dynamicsDynamic wake meandering: calcula�on of wake featuresAeroelas�c wake-aﬀected turbine dynamicsTurbine 1 Turbine 2Meandering wake ﬁeld Figure 4: Schematic of up-wind and wake-a ected wind turbines with a single meandering wake eld. 2.6. Retained sensorial output of aero-servo-elastic simulations Out of the hundreds of sensorial output available from our simulations, we elect to retain only four for the sake of interpretability of the results and brevity of the publication, as shown in Table 4. These include the blade root and the tower base bending moments. Our presumption is to retain the sensorial output of aero-servo-elastic simulations that could, in the real-world, be directly deployed with little technical and economic cost if not already available in today's structural health monitoring systems on WTs in the eld using, for instance, ber-based strain sensing. 9 Table 4: Retained sensorial output of aero-servo-elastic simulations. Sensor Description TwrBsMyt Tower base fore-aft bending moment TwrBsMxt Tower base side-side bending moment RootMyb 1 Blade 1 root apwise bending moment RootMxb 1 Blade 1 root edgewise bending moment 3. Interpreting the simulations output In this section we attempt to provide an interpretation of the numerical simulations, speci cally with respect to the elementary e ects of the in ow RVs on the short-term DEL of the blades and tower base in relation to the wake formation and transport. A qualitative analysis is here o ered with the purpose of informing the subsequent GPR and sensitivity analysis. 3.1. E ect of mean wind speed Uand turbulence U In Figure 5 the up-wind wind turbine is operating below rated power, which corresponds to a high thrust coecient and high induction, thus leading into a signi cant wake de cit (lower wind speed in the wake) and an increase in turbulence. When the spacing is set to 3 5Dthe tower base fore-aft DEL of the wake-a ected wind turbine is higher compared to DEL of the up-wind wind turbine. The high thrust coecient and high induction from the up-wind wind turbine lead in reduction of the mean wind speed and an increase in the turbulence intensity of the transported wake towards the down-wind turbine. A wake-a ected wind turbine experiences increased turbulence on the basis of two main contributions. Firstly, small-scale turbulence due the breakdown of tip vortices and turbulence generated by the shear layer in the edges of the wake and secondly from the meandering of the wake de cit relative to the position of the wake-a ected rotor [24]. This increase in turbulence explains why the short-term tower base fore-aft DEL (welded steel) of the wake-a ected wind turbine is higher compared to DEL of the up-wind wind turbine. This does not hold true for WTG separation between 8 11D, where turbulence and wind speed recovers, nor does it hold true for the blade root ap moment DEL. This implies a more pronounced e ect of mean wind speed on the composite blades versus turbulence for the welded steel tower. For spacing above 8Dthe di erence in DEL for both upwind and wake-a ected wind turbine are marginal for both blades and tower structures. For wind speeds that lie above rated >11m=s and for a spacing of 35D, the upwind turbine tower base fore-aft DEL start to exceed that of the wake-a ected wind turbine ( gures omitted for brevity). This implies that the small scale mixing in the wake dominates the large scale e ects of wake meandering, with the mean wind speed rendered the main driver of DEL. This may be attributed to less ecient transport of the wake meandering, as a result of lower thrust coecient and induction on the rotor at higher wind speeds. It is further noted that for increasing mean wind speed, the gap in the blade root apwise DEL reduces gradually between the upwind and wake a ected wind turbine ( gures omitted for brevity). 3.2. E ect of wind shear, Wind shear induces a periodic higher induction e ect on one part of the rotor, which results in loss of symmetry for the wake de cit. This implies that elevated wind shear reduces the overall eciency of the rotor, while creating a less severe wake in the process. This is particularly true for wind speeds above rated. This e ect is well captured in Figure 6(c). The data is ltered tou2[1525]m=s above rated wind speed, where rotor induction is low to start with. The shear exponent is varied in ranges corresponding to 2[0:050:11], 2[0:130:18] and 2 [0:20:31]. When the shear exponent is low and in the narrow range 2[0:050:11] the blade root apwise bending moment DEL of the wake-a ected wind turbine is shown to exceed the DEL of the up-wind wind turbine. When the shear exponent increases, the up-wind turbine exhibits higher DEL with respect to the wake-a ected wind turbine. For wind speeds below rated, i.e., those corresponding to high thrust and high induction, as shown in Figure 6(a), ineciencies due to wind shear start to appear in the tail of the exceedence probabilities (corresponding to U10m=s) for the blade root apwise bending moment DEL of the wake-a ected WT, with a reversal resulting in higher loads for 2[0:050:11] compared to higher shear exponent ranges. 10 0.5 1 1.5 2 2.5 TwrBsMyt[kNm] 10410-310-210-1100Pe Up-wind Wake-affected(a)Spacing 35D 0 2000 4000 6000 8000 RootMyb1[kNm]10-410-310-210-1100Pe Up-wind Wake-affected (b)Spacing 35D 0.5 1 1.5 2 2.5 TwrBsMyt[kNm] 10410-410-310-210-1100Pe Up-wind Wake-affected (c)Spacing 811D 0 2000 4000 6000 8000 10000 RootMyb1[kNm]10-410-310-210-1100Pe Up-wind Wake-affected (d)Spacing 811D Figure 5: Exceedence probabilities for tower base fore-aft bending moment DEL and blade root apwise bending moment DEL conditional on u2[310]m=s. 0 2000 4000 6000 8000 RootMyb1[kNm]10-310-210-1100Pe Up-wind, [0.05 0.11] Wake-affected Up-wind, [0.13 0.18] Wake-affected Up-wind, [0.2 0.31] Wake-affected (a)U2[310]m=s 4000 5000 6000 7000 8000 9000 RootMyb1[kNm]10-310-210-1100Pe Up-wind, [0.05 0.11] Wake-affected Up-wind, [0.13 0.18] Wake-affected Up-wind, [0.2 0.31] Wake-affected (b)U2[1114]m=s 4000 6000 8000 10000 RootMyb1[kNm]10-310-210-1100Pe Up-wind, [0.05 0.11] Wake-affected Up-wind, [0.13 0.18] Wake-affected Up-wind, [0.2 0.31] Wake-affected (c)U2[1525]m=s Figure 6: Exceedence probabilities for blade root apwise bending moment DEL conditional on increasing shear exponent ranges, and spacing 3 5D. 11 3.3. E ect of horizontal in ow skewness, In Figure 7a the up-wind wind turbine is operating below rated power, resulting in signi cant wake de cit, i.e., lower wind speed in the wake, and increase in turbulence. When the spacing amounts to 35Dand the horizontal in ow skewness is negative 2[101], the tower base side-side bending moment DEL of the wake-a ected WT is higher compared to DEL of the up-wind WT. However, this di erence in Figure 7b vanishes once the horizontal in ow skewness becomes positive 2[1 10]. A similar e ect of the horizontal in ow skewness is also observed on the blade root edgewise bending moment DEL as shown in Figure 7c and 7d. 1000 2000 3000 4000 5000 6000 TwrBsMxt[kNm]10-310-210-1100Pe Up-wind Wake-affected (a) 2[101] 1000 2000 3000 4000 5000 TwrBsMxt[kNm]10-210-1100Pe Up-wind Wake-affected (b) 2[1 10] 5500 6000 6500 7000 RootMxb1[kNm]10-310-210-1100Pe Up-wind Wake-affected (c) 2[101] 5500 6000 6500 7000 7500 RootMxb1[kNm]10-310-210-1100Pe Up-wind Wake-affected (d) 2[1 10] Figure 7: Exceedence probabilities for tower base side-side bending moment DEL and blade root edgewise bending moment DEL, conditional on u2[310]m=s, and spacing 35D. 3.4. Recommendations This work aims to establish a data-driven model to predict the loads in the wake-a ected WT, relying on loads measurements extracted from adjacent WTs, operating on availability of sparse structural measurements across the farm. From a structural health monitoring and life cycle as- sessment point of view, the following guidelines/recommendations can be suggested for in ow RVs: We recommend acquiring wind in ow turbulence data with high accuracy and precision, primarily for wind speeds below rated, when the aim lies in yielding a con dent predic- tor/surrogate of the tower base fore-aft DEL loads of the wake-a ected wind turbine. We recommend acquiring wind in ow shear data with high accuracy and precision, primarily for wind speeds corresponding to maximum thrust (i.e. U10m=s) and above rated wind speed (i.e. U>14m=s), when the aim lies in yielding a con dent predictor/surrogate of blade ap DEL loads of the wake-a ected wind turbine. We recommend acquiring horizontal in ow skewness data with high accuracy and precision, primarily for wind speeds below maximum thrust (i.e. U<10m=s), when the aim lies in yielding a con dent predictor/surrogate of the tower base side-side and blade root edgewise DEL loads of the wake-a ected wind turbine. 12 4. Methodological framework 4.1. Gaussian Processes Consider the function f(xxx)2Rof the input vector xxx2Rn. The function f() is referred to as a Gaussian Process (GP) if its value, when sampled on a nite number of inputs XXX= xxx1xxx2xxxN , follows the multivariate normal distribution N((XXX);KKK(XXX;XXX)), with mean (XXX) = Effff(XXX)gand covariance KKK(XXX;XXX) = E (fff(XXX)(XXX))(fff(XXX)(XXX))> , wherefff(XXX) = f(xxx1)f(xxx2)f(xxxN)T[36, Sec. 2.2]. In turn, the mean and covariance are of the form: (XXX) =2 6664(xxx1) (xxx2) ... (xxxN)3 7775KKK(XXX;XXX) =2 6664k(xxx1;xxx1)k(xxx1;xxx2)k(xxx1;xxxN) k(xxx2;xxx1)k(xxx2;xxx2)k(xxx2;xxxN) ............ k(xxxN;xxx1)k(xxxN;xxx2)k(xxxN;xxxN)3 7775(8) wherek(xxxi;xxxj) is the respective covariance function, a symmetric positive de nite function which measures the similarity between the pair of input RVs xxxiandxxxj. The GP is determined by the mean and covariance functions. Hereafter, the mean function is assumed as zero, while the covariance kernel is selected as the squared exponential, which is de ned as follows [36, pp. 83-84]: k(xxx;xxx0) =2 fexp 1 2nX i=1`2 i(xix0 i)2! (9) where2 f:=k(xxx;xxx) is the function variance, and `2 i,i= 1;:::;n are scaling factors for each one of the input dimensions. The scaling factors determine the smoothness of the function on the respective input dimension: a very large value indicates large di erences between adjacent points, thus leading to non-smooth behavior; a very small value indicates signi cant similarity between remote points, and suggests that there are no signi cant variations on the signal. As will be explained later, these values are adjusted to the observed data, while the obtained scaling factors can be used to determine the in uence of the input dimensions on the function outcomes. 4.2. Gaussian Process regression In a Gaussian Process Regression (GPR), the covariance kernel is utilized to associate the function values observed on a set of training points, XXX= xxx1xxx2xxxN , to a test input vector xxx, and thus provide an estimate of the function value on the test point. To this end, it is rst assumed that a set of noisy values of the function are observed yi,i= 1;:::;N , whereyi=f(xxxi)+wi withwia zero-mean normally and independently distributed process, with variance 2 w. The noisy function values, grouped in the vector yyy:= y1y2yNT, and the function value on the test inputf(xxx) are jointly normally distributed variables, as follows [36, Sec. 2.2]: f(xxx) yyy N0 000N1 ;k(xxx;xxx)kkk(xxx;XXX) kkk(XXX;xxx)KKK(XXX;XXX) +2 wIIIN (10) wherekkk(xxx;XXX) =kkkT(XXX;xxx) = k(xxx;xxx1)k(xxx;xxxN) is the cross-covariance between the test and the training input vectors. Then, using the properties of the multivariate normal distribution, the distribution of the function on the test input conditioned on the the training noisy function valuesyyyis also Gaussian, as follows [36, Sec. 2.2]: p(f(xxx)jyyy;XXX) =Nf(xxx);Q(xxx) (11) with conditional mean f(xxx) and variance Q(xxx), which are calculated as follows: f(xxx) =kkk KKK+2 wIIIN1yyy (12a) Q(xxx) =kkkk KKK+2 wIIIN1kkkT  (12b) and wherek:=k(xxx;xxx),kkk:=kkk(xxx;XXX), andKKK:=KKK(XXX;XXX). 4.3. Bayesian approach for adjustment of the hyperparameters of the GPR The performance of the GPR is de ned by the kernel parameters, comprised by 2 f,`2 ifor i= 1;:::;n , and2 w, which are jointly referred to as the hyperparameters P:=f2 w;2 f;`2 1;;`2 ng. Often, the GPR hyperparameters are optimized via maximization of the marginal likelihood , de ned as follows [36, Sec. 2.3]: lnp(yyyjXXX;P) =1 2yyyT KKK+2 wIIIN1yyy1 2 KKK+2 wIIIN N 2ln 2 (13) 13 This is a non-linear optimization problem, which is typically solved by gradient-based non-linear optimization methods with the help of the partial derivatives of the marginal likelihood with re- spect to each one of the hyperparameters. This optimization results in point estimates of the hyperparameters. Contrariwise, Bayesian methods aim at determining a distribution for the hyperparameters given the available data, based on some original assumptions on the hyperparameter distribution. Therefore, Bayesian methods aim at calculating the hyperparameter posterior distribution [38, pp. 12-13]: p(Pjyyy;XXX) =p(yyyjXXX;P)p(P)p1(yyyjXXX) (14) wherep(P) is the prior hyperparameter distribution, which encapsulates any a-priori knowledge of the hyperparameter distribution, and p(yyyjXXX) comprises the model evidence, de ned as follows: p(yyyjXXX) :=Z p(yyyjXXX;P)p(P) dP (15) where represents the space where Pis de ned. In the case of GPR, an analytical expression for the hyperparameter posterior is not possible, due to the non-linear interaction of the hyperparameters with the likelihood. Instead, Markov Chain Monte Carlo (MCMC) methods can be used to obtain a sample of the hyperparameter posterior [38, Ch. 6-7]. In the analysis presented below, the Metropolis-Hastings algorithm is used for this purpose. Further details on the Metropolis-Hastings sampling method can be found [38, Ch. 7]. 4.4. Sensitivity analysis based on the GPR input scaling factors The input scaling factors `2 iassociated with the squared exponential kernel function de ned in Eq. (9) determine the smoothness of the kernel on the respective input dimension xi. Large positive values of`2 iindicate a very rough behavior of the function on the dimension xi[36, pp. 21-22]. This happens because the correlation between values on xiandxi+ x, with xa small increment, drops very fast. Otherwise, values of `2 iclose to zero indicate that the function is essentially at on the dimension xi[36, pp. 21-22]. In this case, the increment  xrequired to produce a signi cant change in the correlation needs to be very large. This property can be used as a way to evaluate the sensitivity of a function approximated by a GPR to each one of the input dimensions. This is shown below. The squared exponential kernel in Eq. (9) can be factorized as follows: k(xxx;xxx0) =2 fnY i=1exp 1 2`2 i(xix0 i)2 =2 fnY i=1ki(xi;x0 i) (16) ki(xi;x0 i) := exp 1 2`2 i(xix0 i)2 and thus the contribution of input xito the kernel can be decoupled. The value  xis here de ned as the increment in xirequired to decrease by a value the maximum covariance value. More precisely  x2R+is the value such that: ki(xi;xi+ x) =ki(xi;xi) (17) with 0<  1 a small positive number. Applying the de nition of ki(xi;x0 i) and using the fact thatki(xi;xi) = 1, then: exp 1 2`2 ix2  = 1 (18) and then, solving for  x, the following value is obtained: x=`1 ip 2ln(1) (19) The value  xcan be interpreted as the distance required to move along the i-th input to decrease the correlation (similarity) between the function values f(xi) andf(xi+x) by the value . If the value of  xis larger than the range of the data on the i-th input, then the desired change in the correlation is not feasible. With close to zero, this resut would be an indicator of atness of the function in the direction xi. In turn, this indicates that the function f(xxx) is insensitive to the inputxi. Note that although  xcan just be interpreted as the inverse of `i, and both variables hold the same information,  xfacilitates understanding of the GPR sensitivity towards a certain input. 14 4.5. Kullback-Leibler divergence for comparison of GPRs Consider two GPRs represented as Ma:=fyyya;XXXa;PagandMb:=fyyyb;XXXb;Pbg, with di erent training data and hyperparameters. Then, it is necessary to evaluate whether or not the prediction of GPRsMaandMbat a test point xxxis the same. Considering that the GPR prediction at the test point follows a Gaussian distribution, it is possible to use the Kullback-Leibler (KL) divergence to compare if both predictive distributions are the same [9, p. 57]. The KL divergence for a Gaussian distribution takes the form: DKL(xxxjMa;Mb) =1 2Qa(xxx) Qb(xxx)+(fa(xxx)fb(xxx))2 Qb(xxx)+ lnQb(xxx) Qa(xxx)1 (20) where fj(xxx) andQj(xxx), withj=fa;bg, are the GPR predictive mean and variance calculated with Eq. (12) for each corresponding model. The global KL divergence of the predictions obtained with both GPRs can be obtained by integrating over the whole domain XRn, as follows: DKL(Ma;Mb) =Z XDKL(xxxjMa;Mb) dxxx (21) while a marginalized KL divergence for input xi,i= 1;:::;n can be obtained by integrating with respect to the remaining inputs, as follows: DKL(xijMa;Mb) =Z XiDKL(xxxjMa;Mb) dxxxi (22) wherexxxirepresents the input vector after eliminating input xi, andXiis its respective space. Evaluation of the integrals in Eqns. (21) and (22) is not analytically tractable, and instead, numerical approximations are required. The construction of the global and marginalized KL divergences in Eqns. (21) and (22) is based on the assumption that there are no cross-correlations in the predictive distribution, or more precisely, E (f(xxx1)f(xxx1))(f(xxx2)f(xxx2))jyyy;XXX = 0 (23) forxxx16=xxx2. Although this assumption does not comply with the de nition of the GP, it largely simpli es the calculation of the global and marginalized KL divergences. 5. Results 5.1. Prediction and analysis of DELs from local EOPs In this initial analysis, GPR models are built to predict the DELs of a single WT component based on locally measured EOPs. This corresponds to the ideal case when all the WTs are fully instrumented and a complete set of wind eld parameters are available on each wind turbine. The objective of this initial analysis is to determine which input variables mostly a ect the DELs and to determine if there are signi cant di erences between the loads in the up-wind and wake-a ected WTs at di erent structural components. Accordingly, in the present case the input vector xxx2R4is made up by the 10 minute averages of the wind speed U:=x1, turbulence intensity :=x2, shear exponent :=x3, and in ow horizontal skewness := x4of the respective WT. In turn, the output ycorresponds to the DELs calculated from the respective 10 minute loads measured either in the root of one of the blades in the edgewise or apwise directions, or in the tower base in the fore-aft or side-to-side directions. For the construction of the regression, the range of the input variables is normalized within the interval [0;1], while the values of the DELs are scaled down by a factor of 104. This normalization is used to enhance the numerical stability of the models. Individual GPRs are built for the DELs obtained on the up-wind and wake-a ected WTs at di erent spacing con gurations and with di erent W ohler exponents. 5.1.1. Hyperparameter identi cation The Bayesian inference approach based on the Metropolis-Hastings (MH) sampling algorithm described in Section 4.3, is used to estimate of the posterior hyperparameter distribution of the GPR predictor based on the available coupled DWM-FAST simulations. To this end, 180 input- output pairs are randomly selected to calculate the GPR's marginal likelihood (Eq. (13)) within the MH sampling loop, while the remaining ones are used for posterior model validation. To ensure even distribution of the training samples, a sampling approach based on clustering of the complete set of inputs is performed. More precisely, an agglomerative hierarchical clustering tree based on Ward's linkage on the Euclidean distance is applied. Subsequently, 180 clusters are built based on 15 the obtained linkage and a single input and its respective output are randomly extracted from each one of them. Independent log-normal distributions are de ned for each one of the hyperparameters. Similarly, independent log-normal distributions are selected as proposal distributions. The parameters of these distributions are summarized in Table 5. Table 5: Settings of the Metropolis-Hastings algorithm for sampling of the GPR hyperparameter posterior Hyperparameter Prior Proposal Kernel varianceln2 fN 1;102 ln2 fjln(2 f)N ln(2 f)k1;0:4 2 f Kernel scalingln`2 iN 0;102 ln`2 ijln(`2 i)N ln(`2 i);0:4 `2 i;i= 1;:::; 4 Noise varianceln2 wN 1;102 ln2 wjln(2 w)N ln(2 w);0:4 2 w Number of Monte-Carlo samples 104. The symbols ( 2 f), (`2 i), and (2 w)indicate the values of the same quantity drawn in the previous iteration of the MH sampling algorithm. Using the previously described set-up, a total of 10 000 samples are simulated with the MH algorithm. Figure 8 displays the distribution of the GPR hyperparameters based on the sample of the hyperparameter posterior drawn with the MH algorithm, obtained with the data from the up-wind WT tower edgewise and apwise DELs. It is rst noted that although the initial values and prior distributions are similar for all the GPR hyperparameters, the posterior distributions converge to di erent intervals. The distributions for the noise variance, kernel variance and the rst three scale parameters ( `2 1: wind speed; `2 2: turbulence intensity; `2 3: shear exponent) are quite narrow,`2 4: horizontal in ow skewness, has wider distributions in both edgewise and apwise DELs. The latter seems to indicate that the horizontal in ow skewness has a reduced e ect in the DELs. Further analysis based on the obtained GRPs is provided in the sequel. Figure 8: Boxplots displaying the distribution of the GPR hyperparameters sampled with the MH algorithm for the DELs in the blade edgewise and apwise directions. W ohler exponent: 9. 5.1.2. Analysis of DELs based on the obtained GPRs Response surfaces. After estimation of the hyperparameter posterior, it is possible to calculate response surfaces of the DEL for any values in the input space. For instance, 1D slices displaying the DEL as a function of single input variables while the remaining ones are kept xed can be calculated with the help of the obtained GPR. Figs. 9 and 10 display slices of the DEL in the blade {edgewise and apwise direction{ extracted from the GPR with Maximum A Posteriori (MAP) hyperparameter estimates on the up-wind and down-wind WTs. On each frame, the remaining inputs are set to the sample median values, while the spacing between WTs is 11 rotor diameters, and the W ohler exponent is m= 9. In both cases, the GPR predictive mean is similar in both the up-wind and wake-a ected WTs, while a signi cant decrement is found in the DEL of the wake- a ected WT at high wind speeds. On the other hand, the con dence intervals are well con ned around the mid-part of the input range, while the dispersion increases towards the boundaries of the input ranges. The increased dispersion is observed due to the reduced number of points towards the boundaries. Figure 11 shows the GPR predictive mean of the DEL in the blade apwise direction as a function of wind speed and turbulence intensity ( and set to their median values, d= 11, andm= 9). DEL estimates with variance higher than 50 times the minimum predictive variance are censored in the displayed surfaces. As observed in Figure 10, the DEL in the blade apwise direction increases as the wind speed and the turbulence intensity do, with the wind speed making 16 Figure 9: Slices of the DEL predictive mean and 99% con dence intervals (CI) based on a GPR with MAP hyperparameter estimates obtained on the blade in the edgewise direction in the up-wind and down-wind WTs. Other inputs are kept at their training set sample median values indicated on each frame. WT spacing: 11 diameters; W ohler exponent: 9. Figure 10: Slices of the DEL predictive mean and 99% con dence intervals (CI) based on a GPR with MAP hyperparameter estimates obtained on the blade in the apwise direction in the up-wind and down-wind WTs. Other inputs are kept at their training set sample median values indicated on each frame. WT spacing: 11 diameters; W ohler exponent: 9. the maximum e ect. The training data points, displayed as red dots in the surface, coincide with the GPR predictive mean, but also are con ned to the non-censored area, which indicates that the predictive variance is well-bounded around the area spanned by the training points. In practice, DEL estimates in points outside the non-censored area can be deemed as unreliable, while at the same time have a low probability of occurrence, according to the joint distribution of Uand. Error analysis. The predictive performance of the GPR is measured in terms of the Normalized Mean Squared Error (NMSE) calculated on the validation set. The NMSE is de ned as: NMSE =PNval i=1(yif(xxxi))2 PNval i=1y2 i(24) 17 Figure 11: Surface displaying the GPR predictive mean of the DEL in the blade apwise direction as a function of wind speed and turbulence intensity, with the shear exponent and horizontal in ow skewness set on their median values. W ohler exponent: 9. Red dots indicate intput/output data points. whereNvalis the number of samples in the validation set. Figure 12 displays the NMSE obtained on the prediction of the DELs in the blade and tower on the up-wind and wake-a ected WTs, with the wake-a ected WT located at 11 rotor diameters. The error gures in all the cases are quite low, indicating a fair predictability in all cases, while the performance in the up-wind and wake-a ected WTs is similar. In the worst case, corresponding to the tower side-to-side loads, the NMSE is under 3%. On the other hand, the best case, corresponding to the blade apwise loads, the NMSE is about two levels of magnitude lower than that in the tower side-to-side loads. For the blade loads, the NMSE appears to mildly increase with the W ohler exponent. This tendency is followed by both up-wind and down-wind WTs. Figure 13 displays the NMSE as a function of the spacing between WTs for the blade and tower loads in the wake-a ected WT. The overall error performance does not evidence large variations at di erent spacing to the values observed in Figure 12. Moreover, the magnitude order of the NMSE is the same as that found in the previous analysis. 5.1.3. Sensitivity analysis based on the GPR hyperparameters A sensitivity analysis based on the increment method, described in Section 4.4, is performed on the hyperparameter sample obtained from the GPR hyperparameter posterior with the MH sampling algorithm. Figure 14 displays the log-increments obtained for the DELs of di erent components in the up-wind WT as a function of the W ohler exponent for a correlation reduction of = 0:001. The top frame of Figure 14 displays the log-increments obtained in the case of the loads in the blade edgewise direction. According to the results, the wind speed requires the lowest variation to perform a change in the DEL, with the turbulence intensity and shear exponent following in 18 Figure 12: NMSE obtained on the prediction of di erent DELs on the up-wind and wake-a ected WTs based on the GPR with MAP hyperparameters for di erent W ohler exponents. Wake-a ected WT located at 11 rotor diameters. Top left: blade edgewise; top right: blade apwise; bottom left: tower side-to-side; bottom right: tower fore-aft. Figure 13: Median NMSE averaged among W ohler exponents as a function of the spacing between WTs in the wake-a ected WT based on the GPR with MAP hyperparameters. Top left: blade edgewise; top right: blade apwise; bottom left: tower side-to-side; bottom right: tower fore-aft. ranking. On the other hand, the horizontal in ow skewness has a much lesser in uence on the variation of the DELs. Considering that the range of the variables is normalized to the range [0 ;1], then the horizontal in ow skewness turns out to be inconspicuous in the DELs. This result can be contrasted with the response slices displayed in Figure 9 for the wind turbine in the up-wind position. In e ect, the wind speed seems to have the most complex in uence in the DELs, with the turbulence intensity and shear exponent introducing a less complex and almost linear variation, while the horizontal in ow skewness seems to have a very reduced in uence in the DELs. On the other hand, contrasting the increments obtained for di erent W ohler exponents, it appears that this variable makes no di erence on the sensitivity of DELs to input variables. A similar interpretation can be provided for the loads in the blade apwise direction, and in the tower side-to-side and fore-aft directions. In the case of the blade apwise direction and tower fore-aft direction, the DELs are also sensitive to wind speed, turbulence and shear exponent, while 19 insensitive to the horizontal in ow skewness. In contrast, the loads in the tower side-to-side direction are insensitive to both the shear exponent and horizontal in ow skewness. In all cases, the W ohler exponent has no in uence on the sensitivity of the DELs to input variables. Figure 14: Median and 90% CI of the log-increments found for each one of the input variables based on the GPR model hyperparameter sample as a function of the W ohler exponent for di erent components of the up-wind WT. From top to bottom: Blade edgewise loads, blade apwise loads, tower side-to-side loads, tower fore-aft loads. Figure 15 provides a similar analysis, this time on the loads of the wake-a ected WT as a function of the spacing between WTs, with the W ohler exponent set to m= 9 in the case of the blade loads and m= 3 in the case of the tower loads. For the blade edgewise direction loads, the spacing seems to bear no in uence on the sensitivity. However, it now appears that the horizontal in ow skewness has an e ect in the DELs of the wake-a ected WT measured in the blade edgewise direction. In the case of the blade apwise loads, the spacing between WTs seems to in uence the DELs, but not signi cantly. For the tower loads, the sensitivity to the shear exponent and horizontal in ow skewness appears to vary with the spacing between WTs. Particularly, for certain spacing the shear exponent has an e ect on the loads while for others it does not. Finally, the e ect of wind speed and turbulence intensity in the loads appears to be stable in the wake-a ected WT regardless of the spacing to the up-wind WT. 5.1.4. Comparison of up-wind and down-wind DELs The obtained GPRs are also used to determine the di erence of the loads in the up-wind and wake-a ected WTs. Pointwise di erences can be evaluated by directly comparing the DEL predic- tions of the GPRs of the up-wind and down-wind WTs on a test input point xxx, simply as: f(xxx) := fa(xxx)fb(xxx) where fa(xxx) and fb(xxx) are the GPR predictive mean of the respective sensor in the up-wind and down-wind WTs. In this sense, DEL di erences are calculated on a grid of input parameter values 20 Figure 15: Median and 90% CI of the log-increments found for each one of the input variables based on the GPR model hyperparameter sample as a function of the spacing between WTs for di erent components of the wake-a ected WT. From top to bottom: Blade edgewise loads, blade apwise loads, tower side-to-side loads, tower fore-aft loads. to compose hypersurfaces. Likewise, pointwise K-L divergences are calculated as in Eq. (20) on a grid of input parameters values, to construct a surface of pointwise K-L divergences. Figure 16 shows the surfaces obtained on the case of the DELs in the blade edgewise direction for di erent values of wind speed and turbulence intensity, with the remaining input variables set to a xed value (shear exponent: 0.3; horizontal in ow skewness: 2.95 degrees, W ohler exponent: 9, spacing between WTs: 11 rotor diameters). The DEL point-wise di erence shown in the left frame of Figure 16 indicate that the predicted DELs tend to be higher in the up-wind wind turbine on low wind speeds, while lower on higher wind speeds. The point-wise K-L divergences shown in the right frame of Figure 16 demonstrate that higher di erences in DEL predictions are found on low turbulence intensity values on almost the complete range of wind speeds. It is noted that although the point-wise di erence on low wind speeds appears to be signi cant, the increased predictive variance on that region due to a lower number of training points makes the di erence on the mean be overwhelmed by the respective variance, and thus the corresponding K-L divergences evaluated on the same region are not as signi cant. On the contrary, larger K-L divergences are found on regions with a larger population of training samples, which makes these results more reliable. The isolated e ect of a single input variable on the K-L divergence can be obtained by marginal- izing with respect to the remaining input variables. This is done in Figure 17, where marginalized K-L divergences, according to Eq. (22), are calculated for the wind speed, turbulence intensity, shear exponent and horizontal in ow skewness for the DELs in the blade edgewise direction for di erent spacing between WTs. Regarding to the wind speed, the main di erence appears for wind speeds between 5 to 15 m/s, while at high wind speeds the di erence appears to be less pronounced. The DELs obtained on turbulence intensities around 15% also display the largest di erences. Similarly, 21 Figure 16: Point-wise comparison between the DELs calculated in the blade {edgewise direction{ in the up-wind and wake-a ected WTs as a function of wind speed and turbulence intensity, based on the corresponding optimal GPR models. Left frame: Point-wise DEL di erence; right frame: Pointwise K-L divergence. Shear exponent: 0.3; horizontal in ow skewness: 2.95 degrees, W ohler exponent: 9, spacing between WTs: 11 rotor diameters. shear exponents close to 0.5 also introduce large di erence on the DELs. Moreover, the horizontal in ow skewness induces larger changes in the DELs as its absolute value increases, as could be expected. The di erence is more pronounced as the spacing between WTs is reduced, although the highest di erence is found when the spacing is 5 rotor diameters. This tendency is consistent on all the input variables. Figure 17: Marginalized K-L divergences between the DELs of up-wind and wake-a ected WTs based on the optimal GPRs, evaluated for the blade edgewise direction loads and for di erent WT spacings. W ohler exponent: 9. A summary of the global K-L divergence, calculated as in Eq. (21), for the DELs in the blade edgewise and apwise directions and in the tower base in the fore-aft and side-to-side directions for di erent WT spacings and W ohler exponents are provided in Figure 18. The tendency of the DELs on the blade edgewise direction {top left frame, Figure 18{ follows that described for the marginalized K-L divergences in Figure 17, namely, higher di erences are found for short WT spacings. The tendency follows on the di erent W ohler exponents with some variability. On the other hand, in the blade apwise direction, the DELs seem to have higher di erences on shorter spacings but also on higher, while the di erence is not as marked on the intermediate spacing values. The di erence seems to increase as the W ohler exponent increases in value. The DELs in the tower base fore-aft direction appear to have higher di erences on short WT spacings, but this seems to not be the case for the DELs in the side-to-side direction, which appear to have a less predictable relation to the spacing between WTs. 5.2. Prediction of DELs in the wake-a ected WT from remotely measured variables In this second analysis, a more realistic scenario, where a measurements from a fully instru- mented WT (here the up-wind WT) are used to estimate the DELs in a secondary WT (the 22 Figure 18: Global K-L divergences between the DELs of up-wind and wake-a ected WTs based on the GPRs on di erent sensors, as a function of the WT spacing and W ohler exponent. Top-left: blade edgewise direction; top-right: blade apwise direction; bottom-left: tower fore-aft direction; bottom-right: tower side-to-side direction. wake-a ected WT). The objective of this analysis is to assess the feasibility of load prediction on a secondary WT based on measurements from a fully instrumented one. To this end, di erent scenar- ios corresponding to di erent input variable combinations are considered, as summarized in Table 6. Scenario 1 corresponds to a still idealized case, where all the wind features of the up-wind WT are available (speed, turbulence intensity, shear exponent and horizontal in ow skewness). Scenarios 2 and 3 re ect more realistic cases, where a few elementary wind eld properties of the up-wind WT are available (wind speed and turbulence) along with the loads on the blade of the same WT. Finally, scenario 4 corresponds to the case when only loads from the up-wind WT are used to predict loads in the wake-a ected WT. The last three scenarios may conform to the conditions of actual wind parks, where one or a few WTs may be fully instrumented, and the objective is then to predict the loads from fully instrumented WTs to adjacent ones. Scen. Input variables (Measured only in up-wind WT) Uupup up upDELBEWDELBFWDELTFADELTSS 1X X X X 2X X X X 3X X X 4 X X X X Table 6: List of variables measured in the up-wind WT used for prediction of DELs in the wake-a ected WT for each one of the remote prediction scenarios. As in the previous analysis, independent GPRs are built for the DELs measured on the blade root and tower base of the wake-a ected WTs at di erent spacings and with di erent W ohler exponents. Moreover, input and output space normalization is performed as in Sec. 5.1, while hyperparameter optimization is carried out with the Bayesian approach based on MH sampling, whose details are provided in Sec. 5.1.1 and in Table 5. A comparison of the performance, in terms of the NMSE (Eq. (24)), obtained on the prediction of DEL in the wake-a ected WT based on local and up-wind EOP measurements is presented in Figure 19. The performance is presented as a function of the spacing between WTs with the W ohler exponent set to 9 for the blade loads and 3 for the tower loads. The overall predictive ability of the local and remote prediction scenarios is quite similar in all the cases. In the case of the blade DELs shown in Figure 19(a)-(b), the best performance is obtained in scenarios 2 and 3, which include the wind speed and the blade loads of the up-wind WT, while the remaining scenarios seem to follow closely. Otherwise, in the case of the tower loads, the error performance in all the cases is roughly the same. These results seem to indicate that it is feasible to predict the loads in the wake-a ected wind turbine based on measurements obtained from an adjacent WT, although in a practical case it should be necessary to include variables such as the wind direction into the regression model. A global sensitivity analysis of all variables included in the remote prediction scenarios, in terms 23 (a) (b) (c) (d) Figure 19: Comparison of the NMSE as a function of the spacing between WTs obtained on the prediction of DELs on the down-wind WT based on local EOP measurements and remote measurements corresponding to scenarios 1 to 4 in: (a) blade root edgewise direction, (b) blade root apwise direction, (c) tower base fore-aft direction, (d) tower base side-to-side direction. of the log-increments de ned in Sec. 4.4, is performed with the aid of the MC samples of the hyperparameter posterior obtained at each scenario. The resulting median and 90% con dence intervals are displayed in Figure 20 for each one of the sensors in the blade and tower of the wake- a ected WT. The cases where error bars are not visible correspond to variables not used in the respective scenario. Most notably, the distribution of the log-increments appears to coincide among the di erent scenarios. The exception corresponds to the DELs of the blade root apwise direction, where the sensitivity to wind speed and up-wind DELs changes as other variables are introduced in the analysis. In general, it can be observed that the wind speed, turbulence intensity and the DELs in the blade of the up-wind turbine seem to be better predictors of the loads in the wake-a ected WT. 6. Concluding discussions and outlook We demonstrated how a GPR-based model can be built and calibrated to e ectively predict the loads on a wake a ected wind turbine either from local or remote wind eld sensors and/or loads sensors in adjacent WTs in a wind farm. The maximum prediction error, obtained in the tower base DELs in the fore-aft direction, results in less than 4% for the simulations investigated herein. In the case of the blade DELs, the overall prediction error is less than 1%. Although simple, the test con guration allows us to explore sensor selection (scenario) with the aim to minimize the prediction NMSEs according to the spacing between turbines. According to Table 6 scenario 4 tends to yield the lowest NMSE of the predictions for simultaneously the blade-root and tower 24 (a) (b) (c) (d) Figure 20: Sensitivity analysis in the remote sensing scenarios in terms of the log-increments on each one of the input variables included in scenarios 1 to 4. Circles indicate the median value and error bars indicate the 90% con dence intervals obtained from the MC sample of the posterior hyperparameter distribution in the case of: (a) blade root edgewise direction, (b) blade root apwise direction, (c) tower base fore-aft direction, (d) tower base side-to-side direction. base loads on the wake-a ected WT when the spacing is 11 D, but returns the highest prediction error for tower fore-aft when WT spacing is 5 D. It turns out that introducing wind in ow RVs as further input to the GPR model (as in scenarios 2 and 3) tends to reduce the predictions NMSE on the wake-a ected WT when the spacing is 5 Dor less. The sensitivity analysis carried out on the GPR hyperparameters and the K-L divergence reveals that the wind speed and turbulence intensity are the main drivers of the blade root and tower base loads in a WT receiving the freestream in ow as well as in a wake a ected wind turbine (see Figure 14). The wind shear a ects the blade loads to a lesser extent on the upwind wind turbine, while the horizontal in ow skewness has no noticeable e ects on the upwind WT. The wind shear has an increased signi cance on the wake-a ected WT (see Figure 15), while the horizontal in ow skewness marginally a ects the blade DELs in the edgewise direction and tower side-side direction. We stress that the estimations of the GPR model are trustworthy only within the bounds of the variable space used for model calibration; any extrapolation is either not possible or may lead to unpredictable / wrong results. It is therefore important to ensure that the distributions used for the fatigue load diagnostics are not outside the bounds of validity of the trained GPR model. We note that we chose the average wind speed, turbulence, wind shear and wind in ow skewness as the RVs that describe the ambient freestream wind in ow conditions. There may be additional random variables that are equally suitable for the purpose, such as the atmospheric stability that has a complex dependence on the wind speed, turbulence, and wind shear. While combinations of the aforementioned variables usually correspond to certain stability conditions and can be used as a proxy, it may be more convenient to directly include stability as a variable as demonstrated in 25 [23, 27]. Outlook: We purposely initiate from a simulated experiment adopting a simpli ed, yet ed- ucative, con guration of two wind turbines. In next steps, we intend to look into scalability and generalization of GPR with Bayesian learning for larger predictive tasks in arbitrary wind farm layout, possibly harnessing recent progress in scalable GPs, e.g. scalable GPR algorithm, namely, kernel interpolation for scalable structured Gaussian processes (KISS-GP) along with massively scalable GP (MSGP) [29]. In the past the DWM model performance has been validated against numerous measured load data sets, while in the present study, a surrogate model-based approach is validated only against the DWM-FAST model predictions. However, direct validation of the performance of the proposed GPR approach against measured loads data was not currently possible. This is an obvious potential topic for future research, which we are actively working towards. We elucidated how using the DELs one can quantify signi cant di erences in the dynamic response between WTs receiving freesrteam wind and those wake-a ected WTs. These di erences are more apparent in the mid-range of wind speeds and turbulence intensity (wind speed from 5 to 15 m/s and turbulence intensity from 10 to 25%). This result suggests that the fatigue life on a WT subject to freestream in ow can be quite distinct from that of a WT operating in waked condition. Nonetheless, as wind eld parameters are the main drivers of short-term DELs, it is dicult to predict which of the cases will lead to a long-term reduced fatigue life on a given structural component. Nonetheless, regression models such as the GPR presented here can be used as surrogates of a WT in a wind park and perform MC simulations to estimate fatigue life. This would be a subject of study in a future work. Furthermore, in practice, the long term material degradation of various structural components of the up-wind WT di ers from that of the wake- a ected WT. So as time passes, the GPR model predictions will gradually diverge from the real short-term DEL experienced by the wake-a ected wind turbine. Our GPR approach will have to be modi ed in order to capture such multi-valued and time-dependent relationships using some form of overlapping mixture of probabilistic GPR models with recurrent model parameters updates. Although the WT set-up presented in this work is a rather simpli ed approximation of actual layouts in a wind park, it has provided many insights on the random variables in uencing the DELs on a wake-a ected WT. This set-up also allowed for testing of various combinations of wind eld and loads sensors in the up-wind WT as input for predicting the loads on a wake-a ected wind turbine. In an actual wind park, wind direction may also greatly determine the DELs in a wake- a ected WT. Namely, a wake-a ected WT can change from a fully-waked condition to a freestream condition according to the wind direction and the position of adjacent WT. Also, the severity and aggregation of the wake(s) can change according to the number of adjacent WTs from a speci c wind direction. This case will also be considered in a future work. Finally, in section 3 we attempted to provide a qualitative interpretation of the simulations re- sults. Deploying probabilistic deep machine learning algorithms could be one way to automatically and more systematically learn such complex relationships and dependence structures in high dimen- sional spaces. These algorithms however will potentially lack the interpretability of the underlying physical process but would make up for it in inferring functional relationships which are crucial in virtual diagnostics of wake a ected WT from adjacent WT. The consequence is the realization that a physics-informed deep learning approaches might be the way forward. This case is in fact an active area of research of ours. Acknowledgments The authors acknowledge the support of the European Research Council via the ERC Starting Grant WINDMIL (ERC-2015-StG #679843) on the topic of Smart Monitoring, Inspection and Life- Cycle Assessment of Wind Turbines. The authors would further like to acknowledge the support the ERC Proof-of-Concept Grant ERC-2018-PoC WINDMIL RT-DT (funded under H2020-EU.1.1., grant agreement ID 825833). The authors would like to thank OSMOS Hellas SA for the granting of graphic material appearing in Figure 1.