High Spatial Resolution Simulation of Annual Wind Energy Yield Using Near-Surface Wind Speed Time Series

In this paper a methodology is presented that can be used to model the annual wind energy yield (AEYmod) on a high spatial resolution (50 m ˆ 50 m) grid based on long-term (1979–2010) near-surface wind speed (US) time series measured at 58 stations of the German Weather Service (DWD). The study area for which AEYmod is quantified is the German federal state of Baden-Wuerttemberg. Comparability of the wind speed time series was ensured by gap filling, homogenization and detrending. The US values were extrapolated to the height 100 m (U100m,emp) above ground level (AGL) by the Hellman power law. All U100m,emp time series were then converted to empirical cumulative distribution functions (CDFemp). 67 theoretical cumulative distribution functions (CDF) were fitted to all CDFemp and their goodness of fit (GoF) was evaluated. It turned out that the five-parameter Wakeby distribution (WK5) is universally applicable in the study area. Prior to the least squares boosting (LSBoost)-based modeling of WK5 parameters, 92 predictor variables were obtained from: (i) a digital terrain model (DTM), (ii) the European Centre for Medium-Range Weather Forecasts re-analysis (ERA)-Interim reanalysis wind speed data available at the 850 hPa pressure level (U850hPa), and (iii) the Coordination of Information on the Environment (CORINE) Land Cover (CLC) data. On the basis of predictor importance (PI) and the evaluation of model accuracy, the combination of predictor variables that provides the best discrimination between U100m,emp and the modeled wind speed at 100 m AGL (U100m,mod), was identified. Results from relative PI-evaluation demonstrate that the most important predictor variables are relative elevation (Φ) and topographic exposure (τ) in the main wind direction. Since all WK5 parameters are available, any manufacturer power curve can easily be applied to quantify AEYmod.


Introduction
The world's energy supply is facing multiple challenges.The depletion of conventional fuels is unavoidable [1,2], greenhouse gas emissions from the burning of fossil fuels most significantly contributes to global warming [3,4] and the emissions of air pollutants affect human health [3,5].Although nuclear energy production enables the reduction of carbon dioxide (CO 2 ) emissions [6], nuclear power plants bear great short-and long-term risk of accidents [7].In order to reduce and avoid negative impacts of the current use of energy resources on the environment and human health, alternative forms of energy utilization must be found.
Renewable energies provide a clean, environmentally friendly and health-compatible alternative to fossil energies and nuclear energy [2,5].One major renewable energy resource is the kinetic energy contained in the atmosphere, commonly known as wind energy.The potential for wind energy utilization to play a key role in the future global energy mix is enormous.Wind energy could supply more than 40 times [8] the annual global electricity consumption.Consequently, the wind power generation capacity of the world is growing constantly with an average annual rate of about 30% over the last decade [2].In the European Union, directive 2009/28/EC [9] aims to cover 20% of the primary energy demand by renewable energies in 2020, including wind energy.The leading wind power producer in Europe is Germany.Germany aims to supply at least 30% of the energy consumption in 2020 by renewable energies [4].Yet, in some German federal states the utilization of wind energy is still far from being exhausted.For instance, the ministry of Environment, Climate Protection and the Energy Sector of the southwestern German federal state of Baden-Wuerttemberg plans to increase the share of wind energy in the energy mix from ~1% in 2015 to 10% in 2020 [10].In order to achieve this political target, up to 1200 new wind turbines with an average output power between 2.5 MW and 3.0 MW must be installed in a period of only five years [10].
The first step in the onshore assessment of potential wind turbine sites is to quantify the site-specific atmospheric wind energy resource at the wind turbine hub height (~80-100 m) [11].The wind resource is predetermined by the large-scale atmospheric circulation and modified by characteristics of surface roughness [12] and terrain [13].As a result, the local wind resource can vary significantly over short distances [8].In contrast to this, ground-based measurements of long-term wind speed at the landscape level are rare and only available for heights near the surface (10 m above ground level (AGL)).Because of the high spatiotemporal variability of the local wind resource [14,15], the low number of available near-surface wind speed measurement sites alone often limits the detailed assessment of the site-specific wind resource.
To overcome the problem of the low number of wind speed measurements and the strong influence of surface and terrain characteristics on the local wind resource, one option is highly resolved statistical modeling of wind speed at hub height.However, mapping of average wind speed alone is insufficient [16], since not only the central tendencies of wind speed distributions determine the wind resource.Therefore, fitting an appropriate theoretical wind speed distribution to empirical wind speed distributions is crucial [17].Which theoretical distribution fits empirical wind speed distributions best is currently under discussion [18,19].
Due to the limited availability of wind speed measurement sites in Southwest Germany, a region with highly complex topography and mosaic-like land cover pattern, the goals of this study are (i) the quantification of the annual wind energy yield (AEY) on a high spatial resolution grid and (ii) the identification of the most important factors influencing the local wind resource.

Study Area and Wind Speed Measurements
The study area is the German federal state of Baden-Wuerttemberg (Figure 1).The low mountain ranges Black Forest (length ~150 km, width ~30-50 km, highest elevations >1400 m) and Swabian Alb (length ~180 km, width ~35 km, highest elevations >1000 m) are the most complex topographical features with the strongest impact on the wind resource over the study area [20].The top of the Feldberg (1493 m) is the highest elevation in the study area.Approximately 38% (13,700 km 2 ) of the study area is covered with forests [21].More details about land cover and topographical features in the study area are summarized in [20].
The wind speed database used in this study consists of time series of the daily mean wind speed measured from 1 January 1979 to 31 December 2010 at 58 meteorological stations by the German Weather Service (DWD).The height (h S ) of wind speed measurements varies between 3 m AGL (stations Bad Wildbad-Sommerberg, Isny) and 48 m AGL (station Karlsruhe).Data preparation included gap filling, testing for homogeneity and detrending according to [20].
The median wind speed near the surface values ( r U s ) vary in the range 0.5 m/s (station Triberg) to 7.5 m/s (station Feldberg) (Table 1).To extend the database, four measurement stations located in the bordering federal states Hesse and Bavaria were included in this study.Out of the 58 wind speed time series, 48 time series were put into a parameterization dataset (DS1).The remaining 10 time series, for which the original length was less than 10 years, belong to the validation dataset (DS2).

Wind Speed Extrapolation
All wind speed near the surface (U S ) time series were extrapolated to 100 m AGL using the Hellman power law [22][23][24].It was demonstrated by [11] that the power law performs well compared to similar wind speed extrapolation methods.According to [22], the accuracy of the power law increases when stratification effects and the influence of the wind speed are considered.Therefore, the Hellmann exponent (E) was computed on a daily basis.
As has previously been done by [20,25], daily mean wind speed at the 850 hPa pressure level (U 850hPa ) and the height of the 850 hPa pressure level AGL (h 850hPa ), both available from the European Centre for Medium-Range Weather Forecast [26], were used to calculate daily, station-specific E-values: After the E-values were determined, daily, station-specific U S -values were extrapolated to 100 m AGL yielding U 100m,emp :

Probability Distribution Fitting
Prior to the probability distribution fitting, U 100m,emp time series were transformed to empirical cumulative distribution functions (CDF emp ).Afterwards, 67 CDF were fitted to each CDF emp .The goodness of fit (GoF) of each CDF was quantified by calculating the coefficient of determination (R 2 ) from probability plots [19,27] and the Kolmogorov-Smirnov statistic (D) [28][29][30] to the fits.The D-values were obtained by measuring the largest vertical difference between CDF and CDF emp .The transformation of time series, fitting and GoF evaluation were done by EasyFit software (Version 5.5, MathWave Technologies, Dnepropetrovsk, Ukraine) and Matlab ® Software Optimization Toolbox (Release 2015a; The Math Works Inc., Natick, MA, USA).
According to Dand R 2 -value evaluation, which will be presented in detail in the results section, the five-parameter Wakeby distribution (WK5) [31] is clearly the best-fitting distribution.It can be defined by its quantile function [20,25,31,32]: where F is the cumulative probability with U 100m,distr (F) being the associated wind speed value.The four parameters α, β, γ, and δ are distribution parameters and the fifth parameter, ε, is the location parameter.WK5 can be interpreted as a mixed distribution [33] consisting of a left and right part [31,32].This enables WK5 to reproduce shapes of wind speed distributions that other distributions cannot reproduce [25,31].

Predictor Variable Building
A total number of 92 predictor variables (50 m ˆ50 m) covering the study area were built by using the ArcGIS ® 10.2 software (Esri, Redlands, CA, USA).All predictor variables originate from a digital terrain model (DTM), CORINE Land Cover (CLC) data [34] or ERA-Interim reanalysis U 850hPa [26].
The DTM was used to map Φ, τ [35,36], curvature, aspect and slope.The Φ-values were calculated by subtracting the mean elevation of an outer circle around each grid point from the grid point-specific elevation.Five different Φ variants with outer-circle radii of 250 m, 500 m, 1000 m (Φ 1000m ), 2500 m (Φ 2500m ) and 5000 m (Φ 5000m ) were created.
The τ-maps were built for eight main compass directions (northeast (22.5  (247.5 ˝-292.4˝), northwest (292.5 ˝-337.4˝), north(337.5˝-22.4 ˝)) at 200 m radius intervals.This was done by summing angles up to a distance limited to 1000 m.Curvature, aspect and slope were calculated by using the Spatial Analyst Toolbox in ArcGIS.
Roughness length (z 0 ) was derived from CLC data with an original spatial resolution of 100 m ˆ100 m.Roughness length values were assigned to land cover types according to [20] yielding the local roughness length (z 0,l ).Additionally, "effective" roughness length values (z 0,eff ) for the eight main compass directions were calculated.This was done for four different radii around each grid point (100 m, 200 m, 300 m, 400 m).In the end, all z 0 -values were interpolated to 50 m ˆ50 m resolution grids.

Wakeby Parameter Estimation and Modeling
The procedure applied to obtain the Wakeby parameters at every grid point in the study area comprised the following work steps: (1) estimating the Wakeby parameters of every CDF emp based on L-moments [38,39]; (2) analyzing the obtained Wakeby parameters and identifying common characteristics of all distributions ; and (3) modeling target variables (Y) that enable the calculation of all WK5 parameters at every grid point in the study area.To make the WK5 parameter modeling more robust, the WK5 parameters estimated by L-moments were modeled indirectly according to [20,25].
Analyzing the estimated distributions led to the following parameter modeling and calculation approach: First, the estimated left-hand tail of WK5 (Y L ), which is represented by α, β and ε, was modeled: The estimated location parameter ε, which represents the lower bound of the distribution, was directly modeled.Because the L-moment-based WK5 parameter estimation showed that α = 10 at nearly all stations, it was set to this value.The use of a fixed α-value enabled the subsequent calculation of β.
Since Y L affected WK5 parameter estimation up to F = 0.25, exactly as described by [31,32], the percentiles F = 0.30 (Y R1 ), F = 0.50 (Y R2 ), F = 0.75 (Y R3 ) and F = 0.99 (Y R4 ) were modeled to build the right-hand tail of WK5 (Y R ).A system of non-linear equations was solved at every grid point yielding γ and δ: In order to calculate U 100m,mod , Y L and Y R were recombined yielding WK5 with modeled parameters (WK5 mod ).
All Y were computed for every grid point by least squares boosting (LSBoost) [40].This was done by using the Ensemble Learning algorithm LSBoost implemented in the Matlab®Software Statistics Toolbox (Release 2015a; The Math Works Inc.).LSBoost is basically a sequence of simple regression trees, which are called weak learners (B).The objective of LSBoost is to minimize the mean squared error (MSE) between Y and the aggregated prediction of the weak learners (Y pred ).In the beginning, the median of the target variables ( r Y) is calculated.Afterwards, multiple regression trees B 1 , . . ., B m are combined in a weighted manner [41] to improve model accuracy.The individual regression trees are a function of selected predictor variables (X): with p m being the weight for model m, M is the total number of weak learners, and v with 0 < v ď 1 is the learning rate [20,41].
The predictor variable selection process comprised several steps.First, the most appropriate length of outer-circle radii for τ and z 0,eff were determined by the correlation coefficient (r) between τ respectively z 0,eff and Y. Secondly, the importance of the remaining predictor variables was evaluated by predictor importance (PI) which quantifies the relative contribution of individual predictor variables to the model output [21].The PI-values were determined by summing up changes in MSE due to splits on every predictor and dividing the sum by the number of branch nodes.All predictor variables with PI = 0.00 were sorted out.
After PI-evaluation, combinations of predictor variables were tested for their predictive power.Starting with one predictor variable, further predictor variables were added to the model and kept when the model accuracy measures R 2 , mean error (ME), mean absolute error (MAE), MSE and mean absolute percentage error (MAPE) improved [42][43][44].For model parameterization, DS1 data were used.Model validation was done with both DS1 and DS2 data.
Multicollinearity among the predictor variables was investigated by assessing the variance inflation and the condition index in combination with variance decomposition proportions according to [45].

Annual Wind Energy Yield Estimation
The relationship between wind speed and the electrical power output (P) of wind turbines is typically established by a power curve [46].Power curve values are developed from field measurements and can be used for studies involving energy calculations [47].There are three important points characterizing a typical power curve (Figure 2): (1) at the cut-in speed the wind turbine starts to generate usable power; (2) after exceeding the rated output speed the maximum output power (rated power) is generated; and (3) after exceeding the cut-out speed turbines cease power generation and shut down [46].A standard 2.5 MW power curve [48] for onshore wind power plants was applied to calculate the AEY.
Energies 2016, 9, 344 6 of 20 Statistics Toolbox (Release 2015a; The Math Works Inc.).LSBoost is basically a sequence of simple regression trees, which are called weak learners (B).The objective of LSBoost is to minimize the mean squared error (MSE) between Y and the aggregated prediction of the weak learners (Ypred).In the beginning, the median of the target variables ( Y ~) is calculated.Afterwards, multiple regression trees B1, …, Bm are combined in a weighted manner [41] to improve model accuracy.The individual regression trees are a function of selected predictor variables (X): with pm being the weight for model m, M is the total number of weak learners, and v with 0 < v ≤ 1 is the learning rate [20,41].
The predictor variable selection process comprised several steps.First, the most appropriate length of outer-circle radii for τ and z0,eff were determined by the correlation coefficient (r) between τ respectively z0,eff and Y. Secondly, the importance of the remaining predictor variables was evaluated by predictor importance (PI) which quantifies the relative contribution of individual predictor variables to the model output [21].The PI-values were determined by summing up changes in MSE due to splits on every predictor and dividing the sum by the number of branch nodes.All predictor variables with PI = 0.00 were sorted out.
After PI-evaluation, combinations of predictor variables were tested for their predictive power.Starting with one predictor variable, further predictor variables were added to the model and kept when the model accuracy measures R 2 , mean error (ME), mean absolute error (MAE), MSE and mean absolute percentage error (MAPE) improved [42][43][44].For model parameterization, DS1 data were used.Model validation was done with both DS1 and DS2 data.
Multicollinearity among the predictor variables was investigated by assessing the variance inflation and the condition index in combination with variance decomposition proportions according to [45].

Annual Wind Energy Yield Estimation
The relationship between wind speed and the electrical power output (P) of wind turbines is typically established by a power curve [46].Power curve values are developed from field measurements and can be used for studies involving energy calculations [47].There are three important points characterizing a typical power curve (Figure 2): (1) at the cut-in speed the wind turbine starts to generate usable power; (2) after exceeding the rated output speed the maximum output power (rated power) is generated; and (3) after exceeding the cut-out speed turbines cease power generation and shut down [46].A standard 2.5 MW power curve [48] for onshore wind power plants was applied to calculate the AEY.The discrete P-values from the manufacturer power curve were interpolated by a spline to obtain a continuous power curve.The basic attributes of the applied power curve are: cut-in speed U 100m = 3.0 m/s; cut-out speed U 100m = 25.0 m/s; rated output speed U 100m = 13.0 m/s and; rated output power P = 2580 kW.The empirical annual wind energy yield (AEY emp ) was calculated for each station in DS1 and DS2 following [49]: PpU 100m,emp,i q{Z 1 (7) with N = 11,688 being the total number of days in the investigation period and the number of years in the investigation period (Z 1 ).
The average electrical power output (P) was calculated according to [19,50]: The above equation describes the electrical power produced at each wind speed class multiplied by the probability of the specified wind speed class and integrated over all possible wind speed classes [50] with f(U 100m,mod ) being the probability density of U 100m,mod .After P is calculated modeled annual wind energy yield (AEY mod ) can be computed by multiplying P with the respective number of days per year (Z 2 ): AEY mod " P ˆZ2

Summary of the Methodology
The methodology for the quantification of AEY in the study area is summarized in Figure 3.The basic steps are: The discrete P-values from the manufacturer power curve were interpolated by a spline to obtain a continuous power curve.The basic attributes of the applied power curve are: cut-in speed U100m = 3.0 m/s; cut-out speed U100m = 25.0 m/s; rated output speed U100m = 13.0 m/s and; rated output power P = 2580 kW.The empirical annual wind energy yield (AEYemp) was calculated for each station in DS1 and DS2 following [49]: with N = 11,688 being the total number of days in the investigation period and the number of years in the investigation period (Z1).
The average electrical power output ( P ) was calculated according to [19,50]: The above equation describes the electrical power produced at each wind speed class multiplied by the probability of the specified wind speed class and integrated over all possible wind speed classes [50] with f(U100m,mod) being the probability density of U100m,mod.After P is calculated modeled annual wind energy yield (AEYmod) can be computed by multiplying P with the respective number of days per year (Z2):

Summary of the Methodology
The methodology for the quantification of AEY in the study area is summarized in Figure 3.The basic steps are: (1) Extrapolation of near-surface wind speed time series to hub height; (2) Identification of a theoretical distribution that is capable of reproducing various shapes of empirical wind speed distributions; (3) Modeling the estimated parameters of the identified theoretical distribution, based on large-scale airflow, surface roughness and topographic features; (4) Mapping of distribution parameters in the study area; and (5) Calculation of the AEY using a wind turbine-specific power curve.Energies 2016, 9, 344 of 20

Distribution Fitting
According to results from the D-evaluation, WK5 fits 23 CDF emp best.As can be seen in Table 2, the D-value averaged over all stations for WK5 (0.02) is lower than the average D-value of all other theoretical distributions.Another well-fitting distribution is the four-parameter Johnson SB distribution (D = 0.03).The best fitting three-parameter distribution is the inverse Gaussian distribution (D = 0.03).In general, the performance of theoretical distributions defined by three or more parameters is better than the performance of two-and one-parameter distributions.In the case of eight theoretical distributions (Johnson SU, Log-Gamma, Log-Pearson 3, Nakagami, Pareto, Reciprocal, Phased Bi-Exponential, Phased Bi-Weibull) no fit to CDF emp could be achieved and therefore the parameter estimation procedure failed.A widely used theoretical distribution applied to empirical wind speed distributions is the two-parameter Weibull distribution [30,[51][52][53][54][55][56][57][58].However, in this study, the fit of the Weibull distribution is poor (D = 0.10) compared to many other theoretical distributions.These results are in accordance with similar studies where the GoF of various theoretical distributions to empirical Energies 2016, 9, 344 9 of 20 distributions was compared [59][60][61].The Weibull distribution is not even the best-fitting two-parameter distribution, which is the lognormal distribution.The best GoF of a one-parameter distribution was achieved by the also widely used Rayleigh distribution [62,63].However, compared to many distributions defined by more parameters the GoF of the Rayleigh distribution was rather poor (D = 0.10).An explanation for the poor fit of distributions with less than three parameters might be that their capacity for reproducing irregular shapes of empirical distributions is limited.Irregularly shaped empirical wind speed distributions often result from complex topography [64].
The evaluation of averaged R 2 -values confirms results of the D-values evaluation.The best-fitting distribution is WK5 (R 2 = 0.9992), followed by Johnson SB (R 2 = 0.9991).
The superior fit of WK5 is in accordance to GoF measures of empirical near-surface (10 m AGL) wind speed distributions in the study area [20].Based on the results presented in this study it is concluded that WK5 is a universal wind speed distribution for the study area.

Predictor Variable Selection and Importance
The screening of r-values showed that the most appropriate length of outer-circle radius was 1000 m for τ and 200 m for z 0eff .Table 3 lists the predictor variables used for all six least squares boosting models (LSBM) and their relative impact to the model outputs.From the large set of predictor variables, predictor selection finally reduced their number to 14.The main wind directions in the study area are west and southwest.It is therefore reasonable that southwesterly and westerly oriented τand z 0eff -predictor variables have a distinct impact to the model outputs.The highest PI-values for any roughness length predictor variable are found for the LSBM output Y R2 and the western sector (PI = 9.2%).However, the PI-value for Y R1 and the southwestern sector is relatively low (PI = 0.6%).The topographic exposure for the southwestern sector, respectively the western sector, is one of the most important predictor variables for modeling ε, Y L , Y R3 and Y R4 .
It is important to note that U 850hPa was not used to model the left-hand tail of WK5, which represents U 100m,mod -values.Low wind speed values mostly occur when the atmosphere is stably stratified [22].Thus, the influence of U 850hPa on U 100m,mod is rather small.
When modeling Y R3 and Y R4 , the large-scale airflow becomes more important PI = {21.4%,14.8%} because high U 100m,mod -values usually occur when the atmosphere is neutrally stratified [22].
Results from PI-evaluation indicate the fundamental role of relative elevation in wind turbine site assessment.The high PI-values for Φ indicate the great importance of Φ for model outputs.The highest PI-value is 80.5% for Φ 2500m when modeling Y R1 .In contrast, the absolute elevation (ψ) was never used as predictor variable.This is reasonable because sites with high ψ-values are not necessarily exposed to high wind speeds.

Wind Speed Mapping
Median U 100m,mod -values ( r U 100m,mod ) are shown in Figure 4.In large parts (75%) of the study area, r U 100m,mod -values are in the range between 3.0 m/s and 4.0 m/s.In only 0.2% of the study area, r U 100m,mod -values are above 4.9 m/s.Due to the complex topography, high and low r U 100m,mod -values can occur within small distances (<500 m).For example, in the Black Forest, which is characterized by narrow, forested valleys, r U 100m,mod -values are very low.However, there are many exposed mountaintops in close proximity to these valleys where r U 100m,mod -values are high.Beside narrow, forested valleys, lowest r U 100m,mod -values (<3.1 m/s) occur in large cities.In the entire study area the effect of topographic exposure on the modeling results is evident by predominantly higher r U 100m,mod -values at sites exposed to the West and Southwest.

Annual Wind Energy Yield
In Figure 5, the empirical AEY per wind speed class (ΔAEYemp), the modeled AEY per wind speed class (ΔAEYmod), the probability density distributions of WK5mod and the probability density distributions fitted to US-values (US,distr) are presented as a function of wind speed classes (intervals of 0.1 m/s) for the stations Hornisgrinde (Figure 5a) and Laupheim (Figure 5b).
It is clear that percentiles (F = {0.30-0.99})from the right-hand tail of WK5mod contribute more to AEY and are thus more important for the total amount of AEYmod.In Laupheim the mode of U100,mod is 2.3 m/s, whereas highest ΔAEYmod is obtained at 8.0-8.1 m/s.Even at the top of the Hornisgrinde, which is one of the windiest places in the study area, the U100m,mod mode value at 4.2 m/s is clearly lower than the wind speed class assigned to the highest ΔAEYmod-value (9.0-9.1 m/s).Overall, the ΔAEYmod-curves fit ΔAEYemp-values obtained for both stations well.

Annual Wind Energy Yield
In Figure 5, the empirical AEY per wind speed class (∆AEY emp ), the modeled AEY per wind speed class (∆AEY mod ), the probability density distributions of WK5 mod and the probability density distributions fitted to U S -values (U S,distr ) are presented as a function of wind speed classes (intervals of 0.1 m/s) for the stations Hornisgrinde (Figure 5a) and Laupheim (Figure 5b).
It is clear that percentiles (F = {0.30-0.99})from the right-hand tail of WK5 mod contribute more to AEY and are thus more important for the total amount of AEY mod .In Laupheim the mode of U 100,mod is 2.3 m/s, whereas highest ∆AEY mod is obtained at 8.0-8.1 m/s.Even at the top of the Hornisgrinde, which is one of the windiest places in the study area, the U 100m,mod mode value at 4.2 m/s is clearly lower than the wind speed class assigned to the highest ∆AEY mod -value (9.0-9.1 m/s).Overall, the ∆AEY mod -curves fit ∆AEY emp -values obtained for both stations well.The map of AEYmod (Figure 6) shows similar patterns like the 100m,mod U -map.By applying the power curve to U100m,mod, the mean AEYmod-value in the study area is 3.4 GWh/yr.The highest AEYmod-value (13.6 GWh/yr) occurs at the top of the Feldberg.Only in 3% of the study area is AEYmod higher than 5.0 GWh/yr.In 31% of the study area AEYmod is lower than 3.0 GWh/yr with a tendency towards lower AEYmod-values in the southeast, which is mainly due to low U850hPa-values in this part over the study area.In contrast, generally higher AEYmod-values were calculated in the northeast where U850hPa-values are highest at the landscape level.The spatial AEYmod-pattern indicates that the local wind resource is mainly determined by terrain features and surface roughness.The map of AEY mod (Figure 6) shows similar patterns like the r U 100m,mod -map.By applying the power curve to U 100m,mod , the mean AEY mod -value in the study area is 3.4 GWh/yr.The highest AEY mod -value (13.6 GWh/yr) occurs at the top of the Feldberg.Only in 3% of the study area is AEY mod higher than 5.0 GWh/yr.In 31% of the study area AEY mod is lower than 3.0 GWh/yr with a tendency towards lower AEY mod -values in the southeast, which is mainly due to low U 850hPa -values in this part over the study area.In contrast, generally higher AEY mod -values were calculated in the northeast where U 850hPa -values are highest at the landscape level.The spatial AEY mod -pattern indicates that the local wind resource is mainly determined by terrain features and surface roughness.The map of AEYmod (Figure 6) shows similar patterns like the 100m,mod U -map.By applying the power curve to U100m,mod, the mean AEYmod-value in the study area is 3.4 GWh/yr.The highest AEYmod-value (13.6 GWh/yr) occurs at the top of the Feldberg.Only in 3% of the study area is AEYmod higher than 5.0 GWh/yr.In 31% of the study area AEYmod is lower than 3.0 GWh/yr with a tendency towards lower AEYmod-values in the southeast, which is mainly due to low U850hPa-values in this part over the study area.In contrast, generally higher AEYmod-values were calculated in the northeast where U850hPa-values are highest at the landscape level.The spatial AEYmod-pattern indicates that the local wind resource is mainly determined by terrain features and surface roughness.This is underlined by the map extract shown in Figure 7.In the topographically structured Black Forest region, it appears that highest and lowest AEY mod -values occur over horizontal distances shorter than 500 m.This finding is in good accordance to a previous study regarding gust speed in the same area [25].The main wind direction can be inferred from highest AEY mod -values over southwest-facing slopes.This is underlined by the map extract shown in Figure 7.In the topographically structured Black Forest region, it appears that highest and lowest AEYmod-values occur over horizontal distances shorter than 500 m.This finding is in good accordance to a previous study regarding gust speed in the same area [25].The main wind direction can be inferred from highest AEYmod-values over southwest-facing slopes.) for Φ 2500m .The correlation between ψ and AEY mod is relatively weak (r = 0.08).This is due to the fact that some highly elevated Black Forest valleys are sites with the lowest AEY mod -values.The correlation between U 850hPa,0.75 and AEY mod (r = 0.08) is also relatively low since the influence of the large-scale airflow on AEY mod is superimposed by influences of local terrain features and surface roughness.All correlations are highly significant with significance values (p) p ď 0.0000.The exemplary functional relationships between classes of four important predictor variables and AEYmod are shown in Figure 9.The variability of AEYmod-values as a function of U850hPa,0.75(Figure 9a) is lower than the variability of the other displayed predictor variables.This is interpreted to mean that the variability of U850hPa,0.75 is of minor importance for explaining the spatial AEYmod-patterns in the study area.Due to their high roughness, AEYmod is lower over forests and cities (Figure 9b).Areas that are exposed to the southwest (τSW < 2°) show higher AEYmod-values (median: 3.9 GWh/yr) than sheltered areas (τSW > 18°) (median: 1.5 GWh/yr) (Figure 9c).The strongest functional relationship is between Φ2500m and AEYmod (Figure 9d).The assigned median AEYmod-values increase from 1.6 GWh/yr at Φ2500m < −150 m to 4.7 GWh/yr at Φ2500m > 150 m.

Model Validation
The MAPE-values indicate that U100m,mod was simulated accurately (Table 4).They are always below 6% for both DS1 and DS2.The R 2 -values are mostly 0.97 for DS1 percentiles and about 0.95 for DS2 percentiles.The largest downward bias is ME = −0.30m/s for F = 0.99.The exemplary functional relationships between classes of four important predictor variables and AEY mod are shown in Figure 9.The variability of AEY mod -values as a function of U 850hPa,0.75 (Figure 9a) is lower than the variability of the other displayed predictor variables.This is interpreted to mean that the variability of U 850hPa,0.75 is of minor importance for explaining the spatial AEY mod -patterns in the study area.Due to their high roughness, AEY mod is lower over forests and cities (Figure 9b).Areas that are exposed to the southwest (τ SW < 2 ˝) show higher AEY mod -values (median: 3.9 GWh/yr) than sheltered areas (τ SW > 18 ˝) (median: 1.5 GWh/yr) (Figure 9c).The strongest functional relationship is between Φ 2500m and AEY mod (Figure 9d).The assigned median AEY mod -values increase from 1.6 GWh/yr at Φ 2500m < ´150 m to 4.7 GWh/yr at Φ 2500m > 150 m.The exemplary functional relationships between classes of four important predictor variables and AEYmod are shown in Figure 9.The variability of AEYmod-values as a function of U850hPa,0.75(Figure 9a) is lower than the variability of the other displayed predictor variables.This is interpreted to mean that the variability of U850hPa,0.75 is of minor importance for explaining the spatial AEYmod-patterns in the study area.Due to their high roughness, AEYmod is lower over forests and cities (Figure 9b).Areas that are exposed to the southwest (τSW < 2°) show higher AEYmod-values (median: 3.9 GWh/yr) than sheltered areas (τSW > 18°) (median: 1.5 GWh/yr) (Figure 9c).The strongest functional relationship is between Φ2500m and AEYmod (Figure 9d).The assigned median AEYmod-values increase from 1.6 GWh/yr at Φ2500m < −150 m to 4.7 GWh/yr at Φ2500m > 150 m.

Model Validation
The MAPE-values indicate that U100m,mod was simulated accurately (Table 4).They are always below 6% for both DS1 and DS2.The R 2 -values are mostly 0.97 for DS1 percentiles and about 0.95 for DS2 percentiles.The largest downward bias is ME = −0.30m/s for F = 0.99.

Model Validation
The MAPE-values indicate that U 100m,mod was simulated accurately (Table 4).They are always below 6% for both DS1 and DS2.The R 2 -values are mostly 0.97 for DS1 percentiles and about 0.95 for DS2 percentiles.The largest downward bias is ME = ´0.30m/s for F = 0.99.The model performance for DS2 is only marginally worse than for DS1.This indicates the portability of LSBM to other data sets.
Performance measures from the comparison of modeled cumulative distribution functions (CDF mod ) with CDF emp associated with U 100m -time series included in DS2 are shown in Table 5.It appears that the GoF measures for modeled WK5 parameters are better than for many statistical distributions that were directly fitted to CDF emp (compare Table 2).

Conclusions
A methodology is presented that allows assessing the statistical AEY on a high spatial resolution (50 m × 50 m) grid in an area with mosaic-like land cover pattern and complex topography.It was found that highest and lowest AEY occurs in highly textured terrain within very small distances (<500 m).The results of this study therefore emphasize the need to assess AEY at very small spatial scales.This is demonstrated in particular by the great importance of the predictor variables relative elevation and topographic exposure in the main wind direction.
Since the methodology allows for the calculation of all WK5 parameters, the AEY for any manufacturer power curve can be estimated.The methodology is easily portable to other heights above ground level as well as to other study areas.The only requirements for the portability are the availability of the following: (i) near-surface wind speed time series as measured in meteorological networks; (ii) a DTM; (iii) a land cover data set; and (iv) wind speed data not influenced by local topography or land use.
The proposed modeling approach is a useful first step in the exploration of the most appropriate wind turbine sites based on the local wind resource.The produced model outputs and maps are valuable starting points for further in-depth wind turbine site assessment.

Conflicts of Interest:
The author declares no conflict of interest sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Conclusions
A methodology is presented that allows assessing the statistical AEY on a high spatial resolution (50 m ˆ50 m) grid in an area with mosaic-like land cover pattern and complex topography.It was found that highest and lowest AEY occurs in highly textured terrain within very small distances (<500 m).The results of this study therefore emphasize the need to assess AEY at very small spatial scales.This is demonstrated in particular by the great importance of the predictor variables relative elevation and topographic exposure in the main wind direction.
Since the methodology allows for the calculation of all WK5 parameters, the AEY for any manufacturer power curve can be estimated.The methodology is easily portable to other heights above ground level as well as to other study areas.The only requirements for the portability are the availability of the following: (i) near-surface wind speed time series as measured in meteorological networks; (ii) a DTM; (iii) a land cover data set; and (iv) wind speed data not influenced by local topography or land use.
The proposed modeling approach is a useful first step in the exploration of the most appropriate wind turbine sites based on the local wind resource.The produced model outputs and maps are valuable starting points for further in-depth wind turbine site assessment.

Conflicts of Interest:
The author declares no conflict of interest sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Figure 1 .
Figure 1.The study area Baden-Wuerttemberg in Southwest Germany and locations of German Weather Service (DWD) stations.Dots indicate parameterization dataset (DS1) stations; stars indicate validation dataset (DS2) stations.

Figure 1 .
Figure 1.The study area Baden-Wuerttemberg in Southwest Germany and locations of German Weather Service (DWD) stations.Dots indicate parameterization dataset (DS1) stations; stars indicate validation dataset (DS2) stations.

Figure 2 .
Figure 2. Power curve used to calculate empirical annual wind energy yield (AEYemp) and modeled annual wind energy yield (AEYmod) depending on wind speed in 100 m above ground level (AGL) (U100m).Figure 2. Power curve used to calculate empirical annual wind energy yield (AEY emp ) and modeled annual wind energy yield (AEY mod ) depending on wind speed in 100 m above ground level (AGL) (U 100m ).

Figure 2 .
Figure 2. Power curve used to calculate empirical annual wind energy yield (AEYemp) and modeled annual wind energy yield (AEYmod) depending on wind speed in 100 m above ground level (AGL) (U100m).Figure 2. Power curve used to calculate empirical annual wind energy yield (AEY emp ) and modeled annual wind energy yield (AEY mod ) depending on wind speed in 100 m above ground level (AGL) (U 100m ).

( 1 )
Extrapolation of near-surface wind speed time series to hub height; (2) Identification of a theoretical distribution that is capable of reproducing various shapes of empirical wind speed distributions; (3) Modeling the estimated parameters of the identified theoretical distribution, based on large-scale airflow, surface roughness and topographic features; (4) Mapping of distribution parameters in the study area; and (5) Calculation of the AEY using a wind turbine-specific power curve.Energies 2016, 9, 344 7 of 20

Figure 3 .
Figure 3. Schematic representation of the workflow applied to obtain annual wind energy yield (AEY).

Figure 3 .
Figure 3. Schematic representation of the workflow applied to obtain annual wind energy yield (AEY).

) are shown in Figure 4 .
In large parts (75%) of the study area, 100m,mod U -values are in the range between 3.0 m/s and 4.0 m/s.In only 0.2% of the study area, 100m,mod U -values are above 4.9 m/s.Due to the complex topography, high and low 100m,mod U -values can occur within small distances (<500 m).For example, in the Black Forest, which is characterized by narrow, forested valleys, 100m,mod U -values are very low.However, there are many exposed mountaintops in close proximity to these valleys where 100m,mod U -values are high.Beside narrow, forested valleys, lowest 100m,mod U -values (<3.1 m/s) occur in large cities.In the entire study area the effect of topographic exposure on the modeling results is evident by predominantly higher 100m,mod U -values at sites exposed to the West and Southwest.

Figure 4 .
Figure 4. Median of modeled wind speed in 100 m AGL ( 100m,mod U ) in the study area.The legend values indicate highest class values.

Figure 4 .
Figure 4. Median of modeled wind speed in 100 m AGL ( r U 100m,mod ) in the study area.The legend values indicate highest class values.

Figure 5 .
Figure 5. ΔAEYemp and ΔAEYmod as a function of wind speed classes (U) (intervals of 0.1 m/s) as well as the probability of wind speed near the surface fitted to a WK5 distribution (US,distr)-and modeled wind speed in 100 m AGL (U100m,mod)-classes for stations: (a) Hornisgrinde; and (b) Laupheim.

Figure 6 .
Figure 6.AEYmod in the study area.The legend values indicate highest class values.

Figure 5 .
Figure 5. ∆AEY emp and ∆AEY mod as a function of wind speed classes (U) (intervals of 0.1 m/s) as well as the probability of wind speed near the surface fitted to a WK5 distribution (U S,distr )-and modeled wind speed in 100 m AGL (U 100m,mod )-classes for stations: (a) Hornisgrinde; and (b) Laupheim.

Energies 2016, 9 , 344 11 of 20 Figure 5 .
Figure 5. ΔAEYemp and ΔAEYmod as a function of wind speed classes (U) (intervals of 0.1 m/s) as well as the probability of wind speed near the surface fitted to a WK5 distribution (US,distr)-and modeled wind speed in 100 m AGL (U100m,mod)-classes for stations: (a) Hornisgrinde; and (b) Laupheim.

Figure 6 .
Figure 6.AEYmod in the study area.The legend values indicate highest class values.

Figure 6 .
Figure 6.AEY mod in the study area.The legend values indicate highest class values.

Figure 7 .
Figure 7. AEYmod in the Southern Black Forest region.The legend values indicate highest class values.

Figure 8
Figure 8 shows r-values which were calculated between AEYmod and various predictor variables.The r-values confirm the results of the PI-evaluation.The highest and lowest r-values are obtained for the most important predictor variables.The highest absolute r-values are (r = |−0.59|)for τSW and (r = |0.58|)for Φ2500m.The correlation between ψ and AEYmod is relatively weak (r = 0.08).This is due to the fact that some highly elevated Black Forest valleys are sites with the lowest AEYmod-values.The correlation between U850hPa,0.75 and AEYmod (r = 0.08) is also relatively low since the influence of the large-scale airflow on AEYmod is superimposed by influences of local terrain features and surface roughness.All correlations are highly significant with significance values (p) p ≤ 0.0000.

Figure 7 .
Figure 7. AEY mod in the Southern Black Forest region.The legend values indicate highest class values.

Figure 8
Figure 8 shows r-values which were calculated between AEY mod and various predictor variables.The r-values confirm the results of the PI-evaluation.The highest and lowest r-values are obtained for the most important predictor variables.The highest absolute r-values are (r = |´0.59|)for τ SW and (r = |0.58|)for Φ 2500m .The correlation between ψ and AEY mod is relatively weak (r = 0.08).This is due to the fact that some highly elevated Black Forest valleys are sites with the lowest AEY mod -values.The correlation between U 850hPa,0.75 and AEY mod (r = 0.08) is also relatively low since the influence of the large-scale airflow on AEY mod is superimposed by influences of local terrain features and surface roughness.All correlations are highly significant with significance values (p) p ď 0.0000.

Figure 9 .
Figure 9. Boxplots of AEYmod as a function of: (a) 0.75 percentile of the wind speed at the 850 hPa pressure level (U850hPa,0.75);(b) local roughness length (z0,l); (c) topographic exposure in southwest direction (τSW); and (d) relative elevation with outer circle radius of 2500 m (Φ2500m).Boxplot style: red lines indicate medians, boxes indicate interquartile ranges, whiskers indicate 1.5-times interquartile ranges.The legend values indicate highest class values.

Figure 9 .
Figure 9. Boxplots of AEYmod as a function of: (a) 0.75 percentile of the wind speed at the 850 hPa pressure level (U850hPa,0.75);(b) local roughness length (z0,l); (c) topographic exposure in southwest direction (τSW); and (d) relative elevation with outer circle radius of 2500 m (Φ2500m).Boxplot style: red lines indicate medians, boxes indicate interquartile ranges, whiskers indicate 1.5-times interquartile ranges.The legend values indicate highest class values.

Figure 9 .
Figure 9. Boxplots of AEY mod as a function of: (a) 0.75 percentile of the wind speed at the 850 hPa pressure level (U 850hPa,0.75); (b) local roughness length (z 0,l ); (c) topographic exposure in southwest direction (τ SW ); and (d) relative elevation with outer circle radius of 2500 m (Φ 2500m ).Boxplot style: red lines indicate medians, boxes indicate interquartile ranges, whiskers indicate 1.5-times interquartile ranges.The legend values indicate highest class values.

Table 1 .
List of DWD stations and corresponding data features.DS1 stations are indicated by identification numbers (ID) 1-48; DS2 stations are indicated by ID values 101-110.s U ~: median wind speed near the surface values; hS: height.ID Station s U ~ hs ID Station

Table 1 .
List of DWD stations and corresponding data features.DS1 stations are indicated by identification numbers (ID) 1-48; DS2 stations are indicated by ID values 101-110.r U s : median wind speed near the surface values; h S : height.

Table 2 .
Distributions ranked (RK) by Kolmogorov-Smirnov statistic (D)-values with their number of parameters (NP).Dand coefficient of determination (R 2 )-values are averages over all meteorological stations.

Table 3 .
Relative importance of predictor variables used for final least squares boosting models (LSBM) in percent.The top three important predictor variables are highlighted in red.

Table 4 .
Performance measures coefficient of determination (R 2 ), mean error (ME), mean absolute error (MAE), mean squared error (MSE) and mean absolute percentage error (MAPE) calculated from the comparison of empirical and modeled cumulative probabilities (F) associated with U 100m -time series included in DS1 and DS2.

Table 5 .
Performance measures from the comparison of modeled cumulative distribution function (CDF mod ) with empirical cumulative distribution function (CDF emp ) associated with U 100m time series included in DS2.
50dian of modeled wind speed in 100 m AGL U 850hPaWind speed at the 850 hPa pressure level U 850hPa,0.01 1.st percentile of the wind speed at the 850 hPa pressure level U 850hPa,0.3030.thpercentile of the wind speed at the 850 hPa pressure level U 850hPa,0.5050.thpercentile of the wind speed at the 850 hPa pressure level U 850hPa,0.75 75.th percentile of the wind speed at the 850 hPa pressure level U 850hPa,0.99 99.th percentile of the wind speed at the 850 hPa pressure level U SWind speed near the surface