From wind to loads: wind turbine site-specific load estimation using databases with high-fidelity load simulations

We define and demonstrate a procedure for quick assessment of site-specific lifetime fatigue loads, using surrogate models calibrated by means of a database with high-fidelity load simulations. The performance of six surrogate models is assessed by comparing site-specific lifetime fatigue load predictions at ten sites. The surrogate methods include polynomial-chaos expansion, quadratic response surface, universal Kriging, importance sampling, and nearest-neighbor interpolation. Practical bounds for the database and calibration are defined via nine environmental variables, and their relative effects on the fatigue 5 loads are evaluated by means of Sobol sensitivity indices. Of the surrogate-model methods, polynomial-chaos expansion provided an accurate and robust performance in prediction of the different site-specific loads. Although the Kriging approach showed slightly better accuracy, it also demanded more computational resources. Taking into account other useful properties of the polynomial chaos expansion method within the performance comparisons, we consider it to generally be the most useful for quick assessment of site-specific loads. 10 Copyright statement. This article is subject to copyright.


Introduction
Before installing a wind turbine at a particular site, it needs to be ensured that the wind turbine structure is sufficiently robust to withstand the environmentally-induced loads during its entire lifetime.As the design of serially-produced wind turbines is typically based on a specific set of wind conditions, i.e. a site class defined in the IEC (2005) standard, any site where the conditions are more benign than the reference conditions is considered feasible.However, often one or more site-specific parameters will be outside this envelope-and disqualify the site as infeasible, unless it is shown that the design load limits are not going to be violated under site-specific conditions.Such a demonstration requires carrying out simulations over a full design load basis, which adds a significant burden to the site assessment process.
Various methods and procedures have been attempted for simplified load assessment such as expansions based on statistical moments, (Kashef and Winterstein, 1999;Manuel et al., 2001), multivariate regression models of first order (Mouzakis et al., 1999) and second order (Toft et al., 2016), and expansions using orthogonal polynomial basis (Murcia et al., 2016).In the present work, we analyze, refine and expand the existing simplified load assessment methods, and provide a structured approach for practical implementation of a high-fidelity load database for site feasibility assessment.The study aims at fulfilling the following three specific goals: define a simplified load assessment procedure which can take into account all the relevant external parameters required for full characterization of the wind fields used in load simulations; define feasible ranges of variation of the wind-related parameters, dependent on wind turbine rotor size; and obtain estimates of the statistical uncertainty and parameter sensitivities.
The scope of the present study is loads generated under normal power production, which encompasses design load cases (DLC) 1.2 and 1.3 from the IEC 61400-1 standard (IEC, 2005).These load cases are the main contributors to the fatigue limit state (DLC1.2) and often the blade extreme design loads (DLC1.3).The methodology used can easily be applied to other load cases governed by wind conditions with a probabilistic description.Loads generated during fault conditions (e.g.grid drops) or under deterministic wind conditions (e.g.operational gusts without turbulence) will in general not be (wind climate) site-specific.
2 High-fidelity loads database

Definition of variable space
The turbulent wind field serving as input to aeroelastic load simulations can be fully characterized statistically by the following variables: • mean wind field across the rotor plane as described by the -average wind speed at hub height (U ), -vertical wind shear exponent (α), -wind veer (change of mean flow direction with height, ∆ϕ); • turbulence described via -variance of wind fluctuations, σ 2 u , -turbulence probability density function (e.g.Gaussian), -turbulence spectrum defined by the Mann (1994) model with parameters • turbulence length scale L, • anisotropy factor Γ, • turbulence dissipation parameter αε 2/3 ; 2 Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2018-18Manuscript under review for journal Wind Energ.Sci. Discussion started: 2 March 2018 c Author(s) 2018.CC BY 4.0 License.
• air density ρ; • mean wind inflow direction relative to the turbine in terms of -vertical inflow (tilt) angle φv and -horizontal inflow (yaw) angle φh .
All of the parameters above, except the parameters defining mean inflow direction, are probabilistic and site-dependent in nature.The mean inflow direction parameters represent a combination of deterministic factors (i.e.terrain inclination or yaw direction bias in the turbine) and random fluctuations due to e.g.large-scale turbulence structures or variations in atmospheric stability.Mean wind speed, turbulence and wind shear are well known to affect loads and are considered in the IEC61400-1 standard.In Kelly et al. (2014) a conditional relation describing the joint probability of wind speed, turbulence and wind shear was defined.The effect of implementing this wind shear distribution in load simulations was assessed in Dimitrov et al. (2015), showing that wind shear has importance especially for blade deflection.The Mann model parameters L and Γ were also shown to have a noticeable influence on wind turbine loads (Dimitrov et al., 2017).By definition, for given L and Γ the αε 2/3 parameter from the Mann model is directly proportional to σ 2 u L −2/3 (Mann, 1994;Kelly, 2018), and can therefore be omitted from the analysis.The probability density function (pdf) typically used to synthesize time series of velocity components from Mann-model spectra is Gaussian.For a slightly smaller turbine, the NREL 5MW turbine, the assumption of Gaussian turbulence has been shown to not impact the fatigue loads (Berg et al., 2016).The final list of inflow-related parameters thus reads {U, σ u , α, L, ∆ϕ, Γ, φh , φv , ρ}.
The loads experienced by a wind turbine are a function of the wind-derived factors described above, and of the structural properties and control system of the wind turbine.Therefore, a load characterization database taking only wind-related factors into account is going to be turbine-specific.
The variables describing the wind inflow often have a significant correlation between them, and any site-specific load or power assessment has to take this into account using an appropriate description of the joint distribution of input variables.At the same time, most probabilistic models require inputs in terms of a set of independent and identically distributed (i.i.d) variables.
The mapping from the space of i.i.d variables to joint distribution of physical variables requires applying an isoprobabilistic transformation as e.g. the Nataf transform (Liu and Der Kiureghian, 1986), and the Rosenblatt transformation (Rosenblatt, 1952).In the present case, it is most convenient to apply the Rosenblatt transformation, which maps a vector of n dependent variables X into a vector of independent components Y based on conditional relations: . (1) Further mapping of Y to a standard Normal space vector U is sometimes applied, i.e.
For the currently considered set of variables, the Rosenblatt transformation can be applied in the order defined in Table 1 -i.e., the wind speed is considered independent of other variables, the turbulence is dependent on the wind speed, the wind shear is conditional on both wind speed and turbulence, etc.For any variable in the sequence, it is not necessary that it is dependent on all higher-order variables (it may only be conditional on a few of them or even none), but it is required that it is independent from lower-order variables.

Defining the ranges of input variables
The choice for ranges of variation of the input variables needs to ensure a balance between two objectives: a) covering as wide a range of potential sites as possible, while b) ensuring that the load simulations produce valid results.To ensure validity of load simulations, the major assumptions behind the generation of the wind field and computation of aerodynamic forces should not be violated, and the instantaneous wind field should have physically meaningful values.
The turbulence intensity, I u = σ u /U, upper limit can be written as the IEC-prescribed form (ed. 3, sub-class A) with , plus a constant (representing the larger expected range of TI, to span different sites) and a term that encompasses low-windspeed sites and regimes which have higher turbulent intensities.This form is basically equivalent to σ u,IEC + The bounds for turbulence intensity as function of mean wind speed are shown on Figure 2. The limits on shear exponent were chosen following the derivations and findings of Kelly et al. (2014) for P (α|U ), expanding on the established σ α (U ) form to allow for a reasonably wide and inclusive range of expected cases, and also accounting for rotor size per height above ground.This includes an upper bound which allows for enhanced shear due e.g. to lower-level jets and terrain-induced shear; the lower bound also includes the R/z dependence, but does not expand the space to the point that it includes jet-induced negative shear (these are generally found only in the top portion of the rotor).The condition L > max{7.5m,(15m)|α|−2/3 } arises from consideration of the relationship between L, α, σ u , and ε; small shear tends to correlate with larger motions (as in convective well-mixed conditions).The minimum scale (7.5 m) and proportionality constant (15 m) are set to allow a wide range of conditions (though most sites will likely have a scaling factor  larger than 15 m).The maximum Mann-model length scale is chosen based on the limits of where the model can be fitted to measured spectra; this is dictated also by the limits of stationarity in the atmospheric boundary layer.The range of Γ is also dictated by the minimum expected over non-complex terrain within reasonable use of the turbulence model (smaller Γ might occur for spectra fitted at low heights over hills, but such spectra should be modelled in a different way, as in e.g.Mann (2000)).
The range of veer is limited in a way analogous to shear exponent, i.e. it has a basic 1/U dependence; this range also depends upon the rotor size, just as (dU/dz)| rotor = αD/U (Kelly and van der Laan, 2018).The limits for ∆ϕ h above peak follow from the limits on α, while for unstable conditions (∆ϕ h < ∆ϕ h,peak , e.g.all ∆ϕ h < 0) then the veer limit follows a semi-empirical form based on observed extremes of ∂ϕ h /∂z.

Reference high-fidelity load database
Building a large database with high-fidelity load simulations covering the entire variable space is a central task in the present 10 study as such a database can serve several purposes: 1) be directly used as a site assessment tool by probability-weighting the relative contribution of each point to the design loads; 2) serve as an input for calibrating simplified models such as orthogonal-polynomial based expansions; 3) be used as reference for the performance of load models calibrated by other means.
Characterizing the load behaviour of a wind turbine over a range of input conditions requires an experimental design covering the range of variation of all variables with sufficient resolution.In the case of having more than 3-4 dimensions, a full factorial design with multiple levels quickly becomes impractical due to the exponential increase in the number of design points as function of number of dimensions.Therefore, in the present study we resort to a Monte Carlo (MC) simulation as the main approach for covering the joint distribution of wind conditions.For assuring better and faster convergence, we use the lowdiscrepancy Halton sequence in a quasi-Monte Carlo approach (Caflisch, 1998).Figure 1 shows an experimental design based on Halton sequence, compared to crude Monte Carlo and Latin Hypercube designs (Mckay et al., 2000).While a crude Monte Carlo integration has a convergence rate proportional to the square root of the number of samples N , i.e., the mean error ε ∝ N −0.5 , the convergence rate for a pseudo-Monte Carlo with a low-discrepancy sequence results in ε ∝ N −λ , 0.5 ≤ λ ≤ 1.For low number of dimensions and smooth functions, the pseudo-Monte Carlo sequences show a significantly improved performance over the Monte Carlo, e.g.λ → 1, however for multiple dimensions and discontinuous functions the advantage over crude Monte Carlo is reduced (Morokoff and Caflisch, 1995).Nevertheless, even for the full 9-dimensional problem discussed here, it is expected that λ ≈ 0.6 (Morokoff and Caflisch, 1995), which still means about an order of magnitude advantage, e.g., 10 4 pseudo-Monte Carlo samples should result in about the same error as 10 5 crude Monte Carlo samples.

Database specification
A large-scale generic load database is generated using Halton low-discrepancy sample points within the 9-dimensional variable space defined in section 2.3 (Figure 2 shows the bounds for the first 6 variables).The database setup is the following:  Here U is Beta-distributed, while the other variables are uniformly distributed within their ranges.Blue shading shows the site-specific variable distribution for the NKE reference site (site 0, c.f. Table 2/Section 3).
-Up to 10 4 pseudo-random MC sample points are generated, following a low-discrepancy sequence for obtaining evenly distributed points within the parametric space.
-For each sample point, eight simulations, with 1h duration each, are carried out.
-The simulation parameters are tuned to match the required 10-minute statistics (1 h statistics are slightly different due to longer sampling time).
-Each 1h time-series can be split into six 10-minute series, which on average will have the required statistics.
-Simulation conditions are kept stationary over each 1 h simulation period.
-The DTU 10MW reference wind turbine model (Bak et al., 2013), with the basic DTU Wind Energy controller (Hansen and Henriksen, 2013), is used in the HAWC2 aeroelastic software (Larsen and Hansen, 2012).The main quantities of interest from the load simulation output are the short-term fatigue damage-equivalent loads (DEL), and the 10-minute extremes (minimum or maximum, depending on the load type).For each load simulation, four statistics (mean, standard deviation, minimum and maximum values) are calculated for each load channel.For several selected load channels, the 1 Hz DEL for a reference period T ref are estimated using the expression where for obtaining 1 Hz-equivalent DEL over a 1 h period), S i are load range cycles estimated using a rainflow counting algorithm (Rychlik, 1987), and n i are the number of cycles observed in a given range.For a specific material with fatigue properties characterized by an S-N curve of the form K = N • S m (where K is the material-specific Wöhler constant), the fatigue damage D accumulated over one reference period equals

Definition of lifetime damage-equivalent loads
Obtaining site-specific lifetime fatigue loads from a discrete set of simulations requires integrating the short-term damage contributions over the long-term joint distribution of input conditions.The lifetime damage-equivalent fatigue load is defined as where f (X) is the joint distribution of the multidimensional vector of input variables X.With the above definition, S eq,lifetime is a function of the expected value of the short-term equivalent loads conditional on the distribution of environmental variables.
The integration in ( 5) is typically performed numerically over a finite number of realizations drawn from the joint distribution of the input variables, e.g. by setting up a look-up table or carrying out a Monte Carlo simulation.Thus the continuous problem is transformed into a discrete one: where X i , i = 1 . . .N , is the i th realization of X out of N total realizations, and p(X i ) is the relative, discretized probability of X i , which is derived by weighting the joint pdf values of X so that they satisfy the condition Figure 2 shows the distributions of input variables from our high-fidelity database ( §2.4), along with the site-specific distributions for reference site 0 (c.f.Table 2 for site list).
One of the simplest and most straightforward (but not necessarily most precise) ways of carrying out the integrations needed to obtain predicted statistics is to use Importance Sampling ('IS'), where probability weights are applied on each of the database sample points (see e.g.Ditlevsen and Madsen, 1996).The integration (importance sampling) function for determining the expected value of a function g(x) is given by where i = 1 . . .N is the sample point number; ] is a 9-component vector array specifying the values of the 9 environmental variables considered at sample point i; ) is the joint pdf of sample point i according to the site-specific probability distribution; and is the joint pdf of sample point i according to the generic probability distribution used to generate the database for the 9 variables.
Based on the above, it is clear that only points in the database which also have a high probability of occurrence in the sitespecific distribution will have a significant contribution to the lifetime load estimate.Therefore, the IS procedure has relatively slow convergence compared to e.g. a pseudo-MC simulation.Figure 3 shows an example of the convergence of an IS integration for reference site 0, based on using a high-fidelity database with 10 4 points.

Obtaining site-specific results using multi-dimensional interpolation
Estimating an expected function value with a true multi-dimensional interpolation from the high-fidelity database would require finding a set of neighboring points which form a convex polygon.For problem dimensions higher than 3, this is quite challenging due to the non-structured sample distribution.However, it is much easier to find a more crude approximation by simply finding the database point closest to the function evaluation point in a nearest-neighbor approach.This is similar to the table look-up technique often used with structured grids; the denser the distribution of the sample points is, the closer will the results be to an actual Monte Carlo simulation.Finding the nearest neighbor to a function evaluation point requires determining the distances between this point and the rest of the points in the sample space.This is done most consistently in a normalized space, i.e.where the input variables have equal scaling.The cdf (cumulative distribution function) of the variables is an example of such a space, as all cdf's have the same range of (0, 1).and an existing sample is computed as the vector norm of the (e.g.9-dimensional vector) differences between the marginal cdf for the two points: where D = Y − Ŷ is the difference between the current evaluation point Y and the existing sample points in the reference database, Ŷ.The vector ] consists of the marginal cdf functions of the input variables X as obtained using a Rosenblatt transformation.
Since some of the input variables may have significantly bigger influence on the result than other variables, it may be useful to weight the cdf of different variables according to their importance (e.g. by making the weights proportional to the variable sensitivity indices; see Section 4.1.2).

Uncertainty estimation and confidence intervals 10
With the present problem of evaluating the uncertainty in aeroelastic simulations, for any specific combination of environmental conditions, X i , there will be uncertainty in the resulting damage-equivalent loads, S eq (X i ).Part of this uncertainty is statistical by nature and is caused by realization-to-realization variations in the turbulent wind fields used as input to the load simulations.
This uncertainty is normally taken into account by carrying out load simulations with multiple realizations (seeds) of turbulence inputs.Provided that enough load realizations have been generated, the seed-to-seed uncertainty can be characterized by the sample statistics (mean and standard deviation) and with an assumption about the statistical distribution.For the damageequivalent loads which are non-negative by definition, a log-normal distribution is a suitable choice; then the confidence intervals for (S eq,lif etime ) m for confidence level α can be found by: where are the parameters of the log-normal distribution, and Φ denotes the standard Normal cumulative distribution function.It should be noted that the confidence intervals defined above are given in terms of the lifetime DEL raised to the power of the Wøhler slope m, meaning that they actually reflect the range of variation of the lifetime fatigue damage, as visible from eq. ( 4).Using the sample information, confidence intervals can also be determined in a straightforward way using the bootstrapping technique (Efron, 1979).Its main advantage is robustness and no necessity for assuming a statistical distribution of the uncertain variable.
With this approach, each function realization is given an integer index, e.g., from 1 to N for N function realizations.Then, a "bootstrap" sample is created by generating random integers from 1 to N , and, for each random integer, assigning the original sample point with the corresponding index, as part of the new bootstrap sample.Since the generation of random integers allows number repetitions, the bootstrap sample will in most cases differ from the original sample.To obtain a measure of the uncertainty in the original sample, a large number of bootstrap samples are drawn, and the resultant quantity of interest (e.g. the lifetime fatigue load) is computed for each of them.Then, the empirical distribution of the set of outcomes is used to define the confidence intervals.If M bootstrap samples have been drawn, and R is the set of outcomes ranked by value in ascending order, then the bounds for confidence level α equal where the square brackets

Reference sites
The site-specific load calculation methods presented in this study are validated against a set of reference site-specific load calculations on a number of different sites which cover a wide fraction of the variable domain included within the high-fidelity database.In order to show a realistic example of situations where a site-specific load estimation is necessary, the majority of the sites chosen are characterized with conditions which slightly exceed the standard conditions specified by a certain typecertification class.Exceptions are site 0 which has the most measured variables available and is therefore chosen as a main reference site, and the "sites" representing standard IEC class conditions.The IEC classes are included as test sites as they are described by only one independent variable (mean wind speed).They are useful test conditions as it may be challenging to correctly predict loads as function of only one variable using a model based on up to 9 random variables.The list of test sites is given in Table 2.For each site, the joint distributions of all variables are defined in terms of conditional dependencies, and generating simulations of site-specific conditions is carried out using the Rosenblatt transformation, (1).The conditional dependencies are described in terms of functional relationships between the governing variable and the distribution parameters of the dependent variable, e.g. the mean and standard deviation of the turbulence are modelled as linearly dependent on the wind speed as recommended by the IEC 61400-1 standard, while the mean wind shear is dependent on the mean wind speed and on the turbulence, as defined by (Kelly et al., 2014).With this procedure, Pseudo-Monte Carlo samples of the environmental conditions at each site are generated from the respective joint distribution, and fed as input to load simulations.The resulting reference lifetime equivalent loads are then defined as the sample means from the Monte Carlo simulations, while the uncertainty in the lifetime loads is estimated using bootstrapping.

Reduced order models
In this section we present three different reduced order models: 1) Polynomial chaos expansion, 2) Universal Kriging, and 3) Quadratic response surface.All three methodologies are to be calibrated to the database presented in Section 2.4 in order to predict site-specific loads.

Polynomial chaos expansion
Polynomial Chaos Expansion (PCE) is a popular method for approximating a stochastic function of multiple random variables using an orthogonal polynomial basis.In the classical definition of PCE (Ghanem and Spanos, 1991) the input random variables X are defined on (−∞, ∞), with Hermite polynomials typically used as the polynomial basis.1 Choosing a polynomial basis which is orthogonal to a non-Gaussian probability measure turns the PCE problem into the so-called Wiener-Askey or Generalized chaos, (Xiu and Karniadakis, 2002).For the present problem, a Generalized PCE using Legendre polynomials is 10 considered most suitable as the Legendre polynomials P n (ξ) are orthogonal with respect to a uniform probability measure in the interval ξ = [−1, 1], which means that the PCE can conveniently be applied on the cumulative distribution functions of the  variables X which are defined in the interval [0, 1] so that   where F (X i ) is the cumulative distribution function of a variable X i ∈ X, i = 1, . . ., M .With this definition, the PCE represents a model applied to a set of transformed variables which due to the applied transformation are independent and identically distributed ('i.i.d.').Note that eq. ( 11) and the evaluation of the cumulative distribution in general does not ensure independence -so an appropriate transformation should also account for the dependence between variables.In the present case, it is convenient to apply the Rosenblatt transformation as defined in eq. ( 1).For the current implementation of PCE, only eq. ( 1) is required since the expansion is based on the Legendre polynomials, however the transformation to standard Normal space in eq. ( 2) is used for other procedures, e.g. the response surface model discussed later.
Using the notation defined by Sudret (2008), we consider the family of univariate Legendre polynomials P n (ξ).A multivariate, generalized PCE with M dimensions and maximum polynomial degree p is defined as the product of univariate Legendre polynomials where the maximum degree is less than or equal to p.The univariate polynomial family for dimension i can be defined as The multivariate polynomial of dimension M is then defined as The total number of polynomials of this type is The aim of using PCE is to represent a scalar quantity S = g(X) in terms of a truncated sequence S(X) + ε where ε is a zero-mean residual term.With this definition, the multivariate generalized PCE of dimension M and maximum degree p is given by where S j ∈ S = [S 1 , . . ., S N p ] are unknown coefficients which need to be determined, and ξ = [ξ 1 , . . .ξ M ] are functions of X as defined in (11).The most straightforward way of determining S is minimizing the variance of the residual ε using a least-squares regression approach: where N p is the number of polynomial coefficients in the PCE and N e is the number of sampling points in the experimental design.For this purpose, a design experiment has to be set up and the so-called design matrix Ψ needs to be constructed: Under the condition that the residuals are (approximately) Normally-distributed, the solution to equation ( 16) is given by with y = g(x (i) ) being a vector with the outcomes of the functional realizations obtained from the design experiment, where The solution of eq. ( 18) requires that the so-called information matrix (Ψ T Ψ) is well-conditioned, which normally requires that the number of collocation points N e is significantly larger than the number of expansion coefficients N p .Subsequently, the problem grows steeply in size when M and p increase.In such situations, it may be desirable to limit the number of active coefficients by carrying out a Least Absolute Shrinkage and Selection Operator (LASSO) regression (Tibshirani, 1996), which regularizes the regression by penalizing the sum of the absolute value of regression coefficients: where λ is a positive regularization parameter; larger values of λ increase the penalty and reduce the absolute sum of the regression coefficients, while λ = 0 is equivalent to ordinary least-squares regression.

Convergence of PCE
We assess the convergence of PCE by calculating the normalized root-mean-square (NRMS) error between a set of observed quantities (i.e.damage-equivalent loads from simulations) y = g(X (i) ), i = 1 . . .N , and the PCE predictions, ỹ = S(X (i) ), i = 1 . . .N , over the same set of N sample points X (i) : where E[y] is the expected value of the observed variable.not explained by the smooth PCE surface.Therefore, the RMSE error achieved under these conditions can be considered near to the smallest possible for the given PCE order and number of dimensions -and further increase in the number of training points or simulation length will not introduce noticeable improvement.However, it should be noted that the number of training points required for such convergence will differ according to the order and dimension of the PCE, and higher order and more dimensions will require more than the ∼ 3000 points which seem sufficient for a PCE of order 6 and dimension 6, as shown on Figure 10.

Sensitivity indices and model reduction
One useful corollary of the orthogonality in the PCE polynomial basis is that the total variance of the expansion can be expressed as the sum of the contributions from individual terms (Sudret, 2008): Each of the terms in the sum in eq. ( 21) represents the contribution of the variables contained in the respective multivariate polynomials Ψ α,j where j = 0 . . .N p − 1.This property can be used for eliminating polynomials which do not contribute significantly to the variance of the output, thus achieving a sparse, more computationally efficient reduced model.By combining the variance truncation and the LASSO regression technique in eq. ( 19), a model reduction of an order of magnitude or more can be achieved.For example, for a 9-dimensional PCE of order 6, using LASSO regularization parameter λ = 1 and retaining the polynomials which have a total variance contribution of 99.5%, resulted in a reduction of the number of polynomials from 5005 to about 200.
Denoting by F i1,...,is the set of all polynomials dependent on a specific combination of input variables (i 1 , . . ., i s ) (and only on them), the sum of variance contributions over F i1,...,is normalized with the total variance represents the PCE-based Sobol' index with respect to variable set F i1,...,is (Sudret, 2008): Based on eq. ( 22) it is also straightforward to obtain the total Sobol indices for a given variable j by summing all SU i1,...,is where j ∈ (i 1 , . . ., i s ).Note that since each variable appears in multiple cross-terms in the expansion, the contributions of some polynomial coefficients are included multiple times in the total Sobol' indices, meaning that the sum of the total indices will typically exceed 1.
The Sobol indices estimated using the above procedure represent the relative contribution to the model variance from variables following the joint input distribution used to calibrate the PCE.In the present case, this distribution would span the uniform variable space of the high-fidelity database defined in Section 2, and the indices will correspond to the load variation within the entire variable ranges as defined in Table 1.It may be more relevant to compute site-specific Sobol indices.
This can be carried out efficiently by a Monte Carlo simulation on the PCE.For number of dimensions equal to M and for N (pseudo) Monte-Carlo samples the required experimental design represents an N × 2M matrix.This is divided into two N × M matrices, A and B.Then, for each dimension i, i = 1 . . .M , a third matrix AB i is created by taking the i th column of AB i equal to the i th column from B, and all other columns taken from A. The load response function, i.e. the PCE, is then evaluated for all three matrices, resulting in three model estimates: f (A), f (B), and f (AB i ).By repeating this for i = 1 . . .M , simulation-based Sobol' sensitivity indices can be estimated as where j = 1 . . .N is the row index in the design matrices A, B, and AB i (Saltelli et al., 2008).

Universal Kriging with polynomial chaos basis functions
Kriging (Sacks et al., 1989;Santher et al., 2003) is a stochastic interpolation technique which assumes the interpolated variable follows a Gaussian process.A Kriging metamodel is described (Sacks et al., 1989) by where σ 2 is the process variance which is assumed to be constant, and R(w, x, θ) is the correlation between Z(x) and Z(w).The hyperparameters θ define the correlation behavior, in terms of e.g. a correlation length.Given a set of points x = [x 1 , x 2 , . . .x N ] where the true function values y = Y (x) are known, the aim is to obtain a model prediction at a new point, x .Based on Gaussian theory, the known values Y (x) and the Kriging predictor Ŷ (x ) will be jointly Gaussian distributed, and Ŷ (x ) will have the following mean and variance (Santher et al., 2003): Here Ψ is the design matrix collecting the terms constituting the basis functions f (x), Ψ ij = f j (x i ) for i = 1 . . .N and j = 1 . . .P ; r(x ) is the vector of cross-correlations between the prediction point x and the known points x; R is the correlation matrix of the known points, R ij = R(x i , x j , θ) for i, j = 1, . . ., N ; and Using the predictor functions above requires determining the regression coefficients (β), the field variance (σ 2 ), and the correlation hyperparameters (θ).A suitable approach is to find the values of β, σ 2 and θ which maximize the likelihood of y, i.e. minimize the model error in a least-squares sense (Lataniotis et al., 2015): Here the hyperparameters, θ, appear within the correlation matrix R. Having set up the design matrix Ψ, the expansion coefficients β can be determined with the least-squares approach, The field variance is obtained as From ( 28) and ( 29) it follows that β and σ 2 can be expressed as functions of θ.Therefore, calibrating the Kriging model amounts to finding the values of θ which maximize the likelihood.By combining eqns.27-29 this leads to the optimization problem θ = arg min For a problem with M dimensions, we assume that the correlation between sample points can be modelled using an anisotropic separable correlation function ( (Sacks et al., 1989;Lataniotis et al., 2015), which assumes a different correlation parameter for each dimension.The total correlation is expressed as the product of the individual one-dimensional correlation functions, The one-dimensional correlation functions are assumed to follow an exponential relation to the distance h = ( The functional form of the mean field f (x) T β is identical to the generalized PCE defined in eq. ( 18), meaning that the PCE is a possible candidate model for the mean in a Kriging interpolation.We adopt this approach and define the Kriging mean as a PCE with properties as described in section 4.1.
The main practical difference between regression-or expansion-type models such as regular PCE and the Kriging approach

Quadratic response surface
A quadratic-polynomial response surface (RS) method based on Central Composite Design (CCD) is a reduced-order model which, among other applications, is also used for wind turbine load prediction (Toft et al., 2016).The procedure involves different maximum order, it was observed that expansions of order 4 or 5 may not be sufficiently accurate for some response channels.This is illustrated in Figure 12 where the prediction of main shaft torsion loads using order 4 and order 6 PC expansion are compared against other methods, and the order 4 calculation shows a significant bias.Therefore, the PC expansion used for reporting the results in this section is based on the same 6-dimensional variable input used with the quadratic response surface, and has a maximum order of 6.For evaluating confidence intervals from the reduced-order models (Kriging, PCE and quadratic response surface), two reduced-order models of each type are calibrated -one for the mean values, and the other for the standard deviations.This allows computing confidence intervals directly by eq. ( 9), or generating a number of realizations for each sampled combination of input variables, and subsequently computing confidence intervals by bootstrapping (eq.( 10)).In the present, we use the latter approach, because it allows bootstrapping simultaneously a random sample of the input variables and the random seed-to-seed variations within each sample.As a result, the obtained confidence interval reflects the 10 combination of seed-to-seed uncertainty and the uncertainty due to finite number of samples from the distribution of the input variables.This approach also allows consistency with the Importance Sampling and Nearest-Neighbor interpolation techniques, where bootstrapping is also used.Since the lifetime fatigue load is in essence an integrated quantity subject to the law of large numbers, the uncertainty in computations based on a random sample of size N will be proportional to 1/ √ N .Comparing uncertainties and confidence intervals as defined in Section 2.5.5 will therefore only be meaningful when approximately the same number of samples is used for all calculation methods.This approach is used for generating Figures 13 and 14, where the performance of all site-specific load estimation methods is compared for reference site 0, for 8 load channels in total.Figure 13 shows results for tower base and tower top fore-aft and side-side bending moments, and Figure 14 displays the tower top yaw moment, the main shaft torsion, and blade root flapwise and edgewise bending moments.
The results for Site 0 show that for all methods the prediction of blade root and tower top loads is more accurate than the prediction of tower base loads.Also, overall the predictions from the reduced models -the quadratic RS and the PCE, as well as from the Kriging model, are more robust than the IS and nearest-neighbor (NN) interpolation techniques.Similar performance is observed for most other validation sites.The full site-specific results for all load estimation methods are shown in Tables 3,   4, 5, 6, and 7, for the PCE, Kriging, quadratic RS, Importance Sampling, and NN-interpolation techniques respectively.In all tables, the values represent the mean estimates and are normalized to the results obtained with the direct site-specific Monte Carlo simulations.The largest observed errors amount to ≈ 9% with Kriging, ≈ 10% for the PCE, ≈ 10% for the quadratic RS, ≈ 24% for IS, and ∼ 15-17% for NN-interpolation.Noticeably, the low wind speed, high turbulence site 5 seems as the most difficult for prediction-for most load prediction methods this is the site where the largest error is found.All models except the Kriging also show relatively large errors for the IEC class-based sites.That can be attributed to significantly smaller number of samples used for the ICE-based sites (22 samples where only the wind speed is varied in 1m/s steps).As mentioned above, the statistical uncertainty in the estimation of the lifetime DEL will diminish with increasing number of samples.In addition to this effect, as discussed in section 2.5.3, the uncertainty of the IS model can increase when the site-specific distribution has fewer dimensions than the model.It can be expected that this effect is strongest for the IEC class-based sites, as for them only a single variable -the wind speed -is considered stochastic.

One-to-one comparison and mean squared error
Since the prediction of lifetime fatigue loads is the main purpose of the present study, the performance of the load prediction methods with respect to estimating the lifetime DEL is the main criterion for evaluation.However, the lifetime DEL as an integrated quantity will efficiently identify model bias but may not reveal the magnitude of some uncertainties which result in zero-mean error.As an additional means of comparison we calculate the normalized root-mean square (NRMS) error, defined in eq. ( 20) resulting from a point-by-point comparison of load predictions from a reduced-order model against actual reference values.The reference values are the results from the site-specific aeroelastic load simulations for reference site 0. At each sample point, the reference value y i represents the mean DEL from all turbulence seeds simulated with these conditions.The values of the NRMS error for site 0 for Kriging, RS, and PCE-based load predictions are listed in Table 8. Figure 15 shows a one-to-one comparison for a short sequence of sample points.
The RMS error analysis reveals a slightly different picture.In contrast to the lifetime DEL where the Kriging, PCE and RS methods showed very similar results, the RMS error of the quadratic RS is for some channels about twice the RMS error of the other two approaches.

Variable sensitivities
As described earlier in Section 4.1.2,we consider variable sensitivities (i.e. the influence of input variables on the variance of the outcome) in terms of Sobol indices.By definition the Sobol indices will be dependent on the variance of input variables, meaning that for one and the same model, the Sobol indices will be varying under different (site-specific) input variable distributions.Taking this into account, we evaluate the Sobol indices for the two types of joint variable distributions used in this study -1) a site-specific distribution, and 2) the uniform, bounded joint distribution used to generate the database with high-fidelity load simulations.The Sobol indices for the high-fidelity load database variable range are listed in Table 9, while the indices for the site-specific distribution corresponding to reference site 0 are listed in Table 10.The two tables show similar results -the mean wind speed and the turbulence are the most important factors affecting both fatigue and extreme loads.Other two variables which show a smaller but still noticeable influence are the wind shear α, and the Mann model turbulence length 10 scale L. The effect of wind shear is pronounced mainly for blade root loads which can be explained by the rotation of the blades which, if subjected to wind shear, will experience cyclic changes in wind velocity.The effect of Mann model Γ, veer, yaw, tilt, and air density within the chosen variable ranges seems to be minimal, especially for fatigue loads.

Mean extreme loads
Fatigue loads are the main quantity of interest in the present study, however the methods discussed can also be applied to other types of statistics as long as the load simulations carry the necessary information and cover the input variable space.
For example, the conditions specified in the extreme-turbulence load case scenario DLC1.3 in IEC61400-1 class C and B are within the envelope of the high-fidelity database.DLC1.3 requires taking the mean of the extremes from all seeds simulated 5 at a specific condition -meaning that the requirements for the reduced-order models are very similar to the case with fatigue loads.As a demonstration for the feasibility of these calculations, Figure 16 shows prediction of the mean of extremes for the DLC1.3 load case for IEC class IB, carried out by calibrating a PCE and a quadratic RS according to the same procedures as already defined for the fatigue DEL.6 Discussion and conclusions

Discussion
The previous sections of this paper described the procedure for estimating site-specific lifetime damage-equivalent loads using several simplified model techniques which were applied for ten sites.Based on the site-specific lifetime DEL comparisons, three models showed to be viable (robust and sufficiently accurate) choices for quick site-specific load estimation: the PC expansion, Kriging, and the quadratic RS.When estimating lifetime DEL, these methods showed approximately equal levels of uncertainty.However, in the one-to-one comparisons the quadratic RS model showed larger error, especially for sample points corresponding to more extreme combinations of environmental conditions.This is due to the lower order and the relatively small number of calibration points of the quadratic RS, which means that the model accuracy decreases in the sampling space away from the calibration points, especially if there is any extrapolation.This inaccuracy is reflected in the RMS-error from one-to-one comparisons, but is less obvious in the lifetime fatigue load computations which average out errors with zero mean.The universal Kriging model demonstrated the overall lowest uncertainty both in sample-to-sample comparisons and in lifetime DEL computations.This is to be expected since the Kriging employs an already well-performing model (the PCE), and combines it with an interpolation scheme which reduces the uncertainty even further.However, in most cases the observed improvement over a pure PCE is not significant.This indicates that the sources of the remaining uncertainty are outside the models-e.g. the seed-to-seed turbulence variations: the models being calibrated with turbulence realizations different from the ones used to compute the reference site-specific loads.As a result the trend function has the main contribution to the Kriging estimator and the influence of the Gaussian-field interpolation is minimal.A drawback of the Kriging model with respect to the other techniques is the larger computational demands due to the need of computing correlation matrices and the use of the training sample for each new evaluation.
For all site-specific load assessment methods discussed, the estimations are trustworthy only within the bounds of the variable space used for model calibration -extrapolation is either not possible, or may lead to unpredictable results.It is therefore important to ensure that the site-specific distributions used for load assessment are not outside the bounds of validity of the load estimation model.
The variable bounds presented in this paper are based on a certain degree of consideration of atmospheric physics employed in the relationships between wind speed, turbulence, wind shear, wind veer and turbulence length scale.The primary scope is to encompass the ranges of conditions relevant for fatigue load analysis, and the currently suggested variable bounds include all normal-turbulence (NTM) classes.However, in some situations it may be more practical to choose other bound definitionse.g. for consistency with the IEC 61400-1 definitions of the extreme turbulence models, where the currently suggested bounds do not include ETM class A. For the more advanced methods like PCE and Kriging, there is a practical limitation of the number of training points to be used on a single-computer setup.For a PCE the practical limit is mainly subject to memory availability when assembling and inverting the information matrix, and for a PCE of order 6 and with 9 dimensions this limit is on the order of 1-2•10 4 points on a typical desktop computer.For Kriging, the computing time also plays a role: although a similar number of training points could be stored in memory as for the PCE, the computational time is much longer and the practical limit of training points for most applications is less than for the PCE.However, as it was shown in Sections 4.1 and 5, a training sample of 10 4 points or even half of that should be sufficient for most applications in site-specific load prediction.
Considering the overall merits of the load prediction methods analyzed, the PCE provided an accurate and robust performance.The Kriging approach showed slightly better accuracy but at the expense of increased computational demands.Taking this together with the other useful properties of the PCE -e.g. the orthogonality allowing for creating sparse models or doing variance-based sensitivity analysis, we consider the PCE as the most useful method overall.
In addition to the surrogate models presented in this paper, Artificial Neural Networks (ANNs) are interesting alternative candidates.ANN (see Goodfellow et al., 2016) are machine learning models which are gaining large popularity due to their flexibility and history of successful applications in many problems.It is very likely that a deep Neural Network model can provide similar output quality and performance than the methods described in the present study.This is therefore a matter worth of future research.However, the PCE-based models may have a practical advantage over ANNs due to the "whitebox" features such known variance contributions and the possibility of obtaining analytical derivatives which is important for applications to optimization problems.
The results from the site validations showed that for the majority of sites and load channels, the simplified load assessment techniques can predict the site-specific lifetime fatigue loads to within ≈ 5% accuracy.However, it should be noted that this accuracy is relative to full-fidelity load simulations, and not necessarily to the actual site conditions, where additional uncer-  feasibility assessments, which can in a timely fashion aid the decision on whether to discard a given site as unfeasible, or e.g.make additional high-fidelity computations or more measurements of site conditions.The same procedure, but with additional variables (e.g. 3 variables for wake-induced effects as in (Galinos et al., 2016)) may also be useful as objective function or constraint in farm optimization problems.

Conclusions 5
In the present work we defined and demonstrated a procedure for quick assessment of site-specific lifetime fatigue loads using surrogate models calibrated by means of a database with high-fidelity load simulations.The performance of polynomial chaos expansion, quadratic response surface, universal Kriging, importance sampling, and nearest-neighbor interpolation in predicting site-specific lifetime fatigue loads was assessed.Practical bounds of variation were defined for nine environmental variables and their effect on the lifetime fatigue loads was studied.The study led to the following main conclusions: 10 -The variable sensitivity analysis showed that mean wind speed and turbulence are the factors having the highest influence on fatigue and extreme loads.The Mann turbulence length scale and the wind shear exponent were also found to have an influence, and for the wind shear the effect is more pronounced for rotating components such as blades.Within the studied ranges of variation, the Mann turbulence parameter Γ, wind veer, yaw angle, tilt angle, and air density, were found to have small or negligible effect on the loads.
-The best performing models had errors of less than 5% for most sites and load channels, which is in the same order of magnitude as the variations due to realization-to-realization uncertainty.
-A universal Kriging model with a polynomial chaos expansion used as a trend function achieved the most accurate 5 predictions, but also required the longest computing times.
-A polynomial chaos expansion with Legendre basis polynomials was concluded to be the approach with best overall performance.
-The procedures demonstrated in this study are well suited for carrying out quick site feasibility assessments with a particular wind turbine model.

Where-
R is the rotor radius, D the rotor diameter; α ref,LB = 0.15, α ref,U B = 0.22 are reference wind shear exponents at 15m/s wind speed; -Umax = 25m/s is the upper bound of the wind speed;φ is the reference latitude (here chosen as 50 • ).

Figure 1 .
Figure 1.Comparison of several simulation-based experimental design approaches.Examples show random (Monte Carlo and Latin Hypercube) or pseudorandom (Halton sequence) samples of size 100 drawn from a uniform distribution within a unit hypercube.

Figure 3 .
Figure 3. Convergence of an importance sampling (IS) calculation of the blade root moment from the hi-fi database, towards site-specific lifetime fatigue loads for reference site (site 0).
[x] indicate the integer part of x, and R [x] means the value in R with rank equal to [x].Note that the confidence intervals estimated with the above procedure describe only the statistical uncertainty due to the finite number of random samples, and due to seed-to-seed variation.Narrowing them down by e.g.creating a large number of model realizations does not eliminate other model uncertainties as well as uncertainties in the input variables.Wind Energ.Sci.Discuss., https://doi.org/10.5194/wes-2018-18Manuscript under review for journal Wind Energ.Sci. Discussion started: 2 March 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 5 .
Figure 5. Høvsøre site and mast.Permissions by DTU Wind Energy.

Figure 6 .
Figure 6.The fully instrumented mast at Østerild is the blue dot furthest South/Down on the map.Permissions by DTU Wind Energy.

Figure 7 .
Figure 7. Left: Cup anemometer (left) and wind vane (right) at 244 m height, Østerild test site.Permissions by DTU Wind Energy.
Figure10shows the NRMS error for a PCE of order 6 and with 6 dimensions, as function of the number of samples used to train the PCE, and the hours of load simulations used for each sample point.The NRMS error shown on Figure10is calculated based on a set of 500 pseudo-MC points sampled from the joint pdf of reference site 0, and represents the difference in blade root flapwise DEL observed in each of the 500 points vs. the DEL predicted by a PC expansion trained on a selection of points from the high-fidelity database described in Section 2.

Figure 10 Figure 10 .
Figure10illustrates a general tendency that using a few thousand training samples leads to convergence of the values of the PC coefficients, and the remaining uncertainty is due to seed-to-seed variations and due to the order of the PCE being lower than what is required for providing an exact solution at each sample point.Using longer simulations per sample point does not lead to further reduction in the statistical uncertainty due to seed-to-seed variations -with 4000 training samples, the RMSE error for 1h simulation per sample is almost identical to the error with 8h simulation per sample.The explanation for this observation is that the seed-to-seed variation introduces an uncertainty not only between different simulations within the same sampling point, but also between different sampling points.This uncertainty materializes as an additional variance which is

) 20 Wind
Energ.Sci.Discuss., https://doi.org/10.5194/wes-2018-18Manuscript under review for journal Wind Energ.Sci. Discussion started: 2 March 2018 c Author(s) 2018.CC BY 4.0 License.where x represents the input variables, and Y (x) is the output.The term β T f (x) is the mean value of the Gaussian process (a.k.a. the "trend") represented as a set of basis functions f (x) = [f 1 (x), . . ., f P (x)] and regression coefficients β = [β 1 , . . ., β P ]; Z(x, w) is a stationary, zero-mean Gaussian process.The probability distribution of the Gaussian process is characterized by its covariance
is in the way the training sample is used in the model: in the pure regression-based approaches the training sample is used to only calibrate the regression coefficients, while in Kriging as in other interpolation techniques the training sample is retained and used in every new model evaluation.As a result the Kriging model may have an advantage in accuracy since the model error tends to zero in the vicinity of the training points; however this comes at the expense of an increase in the computational demands for new model evaluations.The extra computational burden is mainly the time necessary to assemble r(x ), the matrix of cross-correlations between the prediction points and the training sample, and the time to multiply r(x ) with other equation terms.Thus, while for a PCE the model evaluation time t(N ) for a sample of size N would follow t(N ) = O(N ), for a Kriging model t(N ) = O(N 2 ).For a Kriging model, a gain in accuracy over the model represented by the trend function willonly materialize in problems where there is a noticeable correlation between the residual values at different training points.In a situation where the model error is independent from point to point (as e.g. in the case when the error is only due to seed-to-seed variations in turbulence) the inferred correlation length will tend to zero and the Kriging estimator will be represented by the trend function alone.

Figure 11 .
Figure 11.Example of a rotatable Central Composite Design (CCD) in a 2-dimensional standard Normal space [u1, u2].The CCD consists of a central point, a 2 k factorial design with 2 levels and k = 2 dimensions, and axial points at distance u = √ 2, meaning that all the outer points lie on a circle.

O6Figure 12 .
Figure 12.Comparison of predictions of the lifetime damage-equivalent loads for six different estimation approaches.All values are normalized with respect to the mean estimate from a site-specific Monte Carlo simulation, and the error bars represent the bounds of the 95% confidence intervals.Results from two PCEs are shown: the blue bar corresponds to the output of a 4 th order PCE, while the black bar corresponds to a 6 th order PCE.
Figure 13.Comparison of predictions of the lifetime damage-equivalent loads for six different estimation approaches.All values are normalized with respect to the mean estimate from a site-specific Monte Carlo simulation.
Figure 14.Comparison of predictions of the lifetime damage-equivalent loads for six different estimation approaches.All values are normalized with respect to the mean estimate from a site-specific Monte Carlo simulation.

Figure 15 .
Figure 15.Sample-by-sample comparison of the 1Hz damage-equivalent load predictions, for two load prediction methods -quadratic Response Surface, and Polynomial Chaos expansion, compared against site-specific Monte Carlo simulation.Each point on the x-axis represents a sample point from the site-specific distribution of input variables for reference site 0, and the y-axis represents the estimated mean 1Hz-equivalent load.

Figure 16 .
Figure 16.Comparison of predictions of the mean of extremes obtained for the extreme turbulence load case DLC1.3 for IEC turbulence class B conditions.MC simulation (blue lines with crosses) denotes the reference values obtained with full-fidelity aeroelastic simulations.

Table 1 .
Bounds of variation of the variables considered.All values are defined as statistics over 10-minute reference period.

Table 2 .
Reference sites used for validation of the site-specific load estimation methods.

Table 3 .
Lifetime damage-equivalent fatigue load predictions obtained using 6-dimensional Polynomial Chaos Expansion of order 6, relative to predictions from site-specific Monte Carlo simulations.Load channel abbreviations are the following: TB: tower base; TT: tower top; Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz:

Table 4 .
Lifetime damage-equivalent fatigue load predictions obtained using the Universal Kriging approach with a 6-dimensional Polynomial Chaos Expansion of order 6 used as basis functions.The results are compared against site-specific Monte Carlo simulations.Load channel abbreviations are the following: TB: tower base; TT: tower top; MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz: torsion.Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My

Table 5 .
Lifetime damage-equivalent fatigue load predictions obtained using 6-dimensional Quadratic Response Surface, relative to predictions from site-specific Monte Carlo simulations.Load channel abbreviations are the following: TB: tower base; TT: tower top; MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz: torsion.Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My

Table 6 .
Lifetime damage-equivalent fatigue load predictions obtained using Importance Sampling of points from a 6-dimensional highfidelity load database, relative to predictions from site-specific Monte Carlo simulations.Load channel abbreviations are the following: TB: tower base; TT: tower top; MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz: torsion.Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My

Table 7 .
Lifetime damage-equivalent fatigue load predictions obtained using nearest-neighbor interpolation in a 6-dimensional high-fidelity load database, relative to predictions from site-specific Monte Carlo simulations.Load channel abbreviations are the following: TB: tower Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My base; TT: tower top; MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz: torsion.

Table 8 .
Root mean square error characterizing the difference between aeroelastic simulations and reduced-order models.Load channel Mx TB My TT Mx TT My TT Mz MS Mz BR Mx BR My abbreviations are the following: TB: tower base; TT: tower top; MS: main shaft; BR: blade root.Loading directions consist of Mx: fore-aft (flapwise) bending, My: side-side (edgewise) bending, and Mz: torsion.

Table 9 .
PCE-based Sobol sensitivity indices for the high-fidelity load database variable ranges.

Table 10 .
Site-specific Sobol sensitivity indices derived for Site 0 using Monte Carlo simulation from a PCE.Tower base fore-aft moment Mx 0.68 0.44 0.03 0.02 0.01 Tower base side-side moment My 0.52 0.40 0.07 0.22 0.17