Estimating vertical wind power density using tower observation and empirical models over varied desert steppe terrain in northern China

A complex and varied terrain has a great impact on the distribution of wind energy resources, resulting in uncertainty in accurately assessing wind energy resources. In this study, three wind speed distributions of kernel, Weibull, and Rayleigh type for estimating average wind power density were first compared by using meteorological tower data from 2018 to 2020 under varied desert steppe terrain contexts in northern China. Then three key parameters of scale factor (c) and shape factor (k) from the Weibull model and surface roughness (z₀) were investigated for estimating wind energy resource. The results show that the Weibull distribution is the most suitable wind speed distribution over that terrain. The scale factor (c) in the Weibull distribution model increases with an increase in height, exhibiting an obvious form of power function, while there were two different forms for the relationship between the shape factor (k) and height: i.e., the reciprocal of the quadratic function and the logarithmic function, respectively. The estimated roughness length (z₀) varied with the withering period, the growing period, and the lush period, which can be represented by the estimated median value in each period. The maximum and minimum values of surface roughness length over the whole period are 0.15 and 0.12 m, respectively. The power-law model and the logarithmic model are used to estimate the average power density values at six specific heights, which show greater differences in autumn and winter, and smaller differences in spring and summer. The gradient of the increase in average power density values with height is largest in autumn and winter, and smallest in spring and summer. Our findings suggest that dynamic changes in three key parameters (c, k, and z₀) should be accurately considered for estimating wind energy resources under varied desert steppe terrain contexts.

1 Introduction

Wind energy is a renewable, environmentally friendly, and popular alternative source of clean energy (Islam et al., 2013; Gabbasa et al., 2013), and as a source of power it has great potential (Chaurasiya et al., 2019). In 2020, 93 GW of new installations brought the global cumulative wind power capacity up to 743 GW. In the onshore market, 86.9 GW was installed, an increase of 59 % compared to 2019. China and the United States remain the world's largest markets for new onshore installations (Joyce and Feng, 2021). To use this kind of nonpolluting energy, a lot of research has been conducted through a variety of different methods to develop an accurate and reliable wind energy evaluation model.

The probability density function (PDF) of wind speed can effectively characterize wind speed. Therefore, the wind speed PDF is of great significance in selecting wind turbine sites, in wind farm design, in generator design, in determining the dominant wind direction, and in evaluating the management and operation of wind conversion systems (Masseran, 2015; Li and Shi, 2010). Wind displays large differences with various topographies, landforms, and meteorological conditions. The magnitude and direction of the wind speed exhibit significant differences when wind flows over rough ground or obstacles in a complex terrain. In addition, the surface topography and roughness of the area around the location of the wind measurement tower will affect the predicted wind resources (Kim and Lim, 2017). Therefore, the wind speed PDF and roughness are important input factors in the estimation of the power density of wind energy.

Different distribution functions have different fitting effects on the actual wind speed values in distinct study areas. According to previous studies (Lo Brano et al., 2011; Celik, 2004; Masseran et al., 2012), seven wind PDFs have been widely used to fit the actual wind speed values: i.e., Weibull, Rayleigh, lognormal, gamma, inverse Gaussian, Pearson type V, and Burr. These models exhibited different advantages and disadvantages for estimating wind probability density. For instance, Celik (2004) used the Weibull and Rayleigh models to perform a statistical analysis of wind energy density in southern Turkey and found that the Weibull model not only fits the measured monthly probability density distribution better but also provides a better power density estimation compared to the Rayleigh model. Masseran et al. (2012) used nine different wind speed PDF models to describe wind speed conditions in different regions of Malaysia and found that gamma, Weibull, and inverse gamma models fit the wind speed data better. Chang (2011a) used six different PDFs, namely, Weibull, a mixture of gamma and Weibull, a mixture of normal, a mixture of normal and Weibull, a mixture of Weibull, and principle of maximum entropy distributions. They were tested on the wind data of three wind farms in Taiwan and it was found that, when the current wind speed distribution is unimodal, the fitting effects of these six PDFs are not significantly different. When the wind speed distribution is bimodal, the other five PDFs are better than Weibull at describing wind characteristics. In addition, many other PDFs have been invented to provide more accurate results for the estimation of wind power density in a specific area (Masseran, 2015; Carta et al., 2009; Jaramillo and Borja, 2004).

Among the aforementioned types of wind speed PDFs, the Weibull and Rayleigh distributions remain the more traditional and widely applicable typical wind speed distribution forms. The key issue in the study of the Weibull distribution is how to accurately determine the values of Weibull scale factor c and shape factor k (Azad et al., 2014; Kaplan, 2017). Generally, six different methods, i.e., graphical method (Basu et al., 2009), empirical method (Costa Rocha et al., 2012; Kaoga et al., 2014), maximum likelihood method (Andrade et al., 2014; Azad et al., 2014), energy trend method (Chang, 2011b; Akdağ and Dinler, 2009), energy pattern method (Andrade et al., 2014), and the moment method (Azad et al., 2014; Kaplan, 2017; Costa Rocha et al., 2012), have been employed to calculate the c and k of the Weibull distribution model. But these methods perform differently in different regions. For instance, Kaplan (2017) found that the energy pattern method and the moment method were the best methods between 2009 and 2013 in the Hatay and Osmaniye regions. When the time series of wind data is provided, the maximum likelihood method is more robust and accurate than other methods (Seguro and Lambert, 2000; George, 2014). In addition, there is a strong time dependence and a high change dependence for the changes in shape factor k and scale factor c (Lun and Lam, 2000; Justus and Mikhail, 1976): e.g., the scale factor c has a power-law functional relationship with height and the shape factor k has a reciprocal logarithmic functional relationship with height. Therefore, we can explore its general laws by studying the seasonal changes and height changes in shape and scale parameters in a specific area.

Roughness length plays a key role in estimating wind energy resources. For example, Laporte (2010) pointed out that the roughness estimation error can cause 5–10 % of the wind energy resource estimation error. Current wind energy resource assessment is based on measured wind data at a height of 60–80 m from the ground, but the actual height of the hub may be greater than these heights. Therefore, we need to combine the surface roughness length and the known wind speed value of the measured height to extrapolate the wind speed value at the height of the hub (Nayyar and Ali, 2020). Theoretically, the surface roughness length z₀ is the height at which the average wind decreases to zero with height. The value z₀ varies with the underlying surface (Davenport et al., 2000; Duan et al., 2021). Currently, three approaches (the analysis method, the Charnock method, and the statistical method) are widely applied to estimate the surface roughness length of offshore wind energy (Golbazi and Archer, 2019). Among them, the statistical method is convenient, as it needs only three layers of wind speed data. After comparing the average value and median value of roughness z₀, it is found that the median value is an order of magnitude closer to the roughness length calculated from the other two methods. Therefore, when using the field measurement method to statistically determine the surface roughness length, attention should be paid to using the median value instead of the average value; otherwise huge errors will be generated when the wind speed is extrapolated to the height of the hub, which will have a major impact on the evaluation of wind energy resources.

As an important production base of wind power energy in northern China, Inner Mongolia is under the influence of the westerly wind all year round. The types of underlying surfaces of wind power towers in China are complex and diverse, including offshore, mountainous, urban outskirts, and grasslands. In Inner Mongolia, especially the desert grassland, the terrain is open, the vegetation is low and sparse, and its wind resources are very rich. Thus, taking the Ningyuanbailiutu site as an example, in-depth data mining was carried out on the four heights of 10, 30, 50, and 70 m for the meteorological element data of a 100 m wind tower from the autumn of 2018 to the summer of 2020 in Damaoqi, Baotou City, Inner Mongolia, China. The following three steps are used to study the three important key parameters that affect the evaluation of wind energy resources: the surface roughness length z₀, the scale factor c, and the shape factor k in the Weibull distribution function. First, we need to determine the uniqueness and importance of the Weibull distribution function in the wind speed time series data in the Damaoqi area. This is reflected in the advances and shortcomings of the kernel distribution model, the Rayleigh distribution model, the Weibull distribution function, and the frequency distribution model using actual wind speed, which are used to calculate the monthly, seasonal, and all-time average power densities. Second, by studying the monthly and seasonal changes in the surface roughness length and the changes in different incoming flow directions, we will gain a comprehensive understanding of the roughness of the site area in Inner Mongolia. Finally, by using two different models, namely, the power-law model with scale parameter c and the logarithmic model with roughness information, the average wind power densities at six specific heights (75, 80, 85, 90, 95, and 100 m) per month, per season, and throughout the period are calculated. In this way, we discuss the application significance of the two models for wind energy development, and provide a scientific reference to further our understanding of the wind energy resources in the region.

2 Study site, data, and methods

2.1 Study site and data

In this study, long-term in situ measurement was conducted in Damaoqi, Baotou City, Inner Mongolia (42^∘04^′25.738^′′ N, 110^∘29^′2.778^′′ E; 1376 m above sea level) from 1 September 2018 to 31 August 2020 (Fig. 1). The observation wind tower is located at the northern foot of Daqing Mountain in the central area of the Inner Mongolia Plateau. Wind speed and wind direction (010C cup anemometers and 020C wind vanes, Metone, USA), atmospheric pressure (CS106 Campbell, USA), air temperature, and humidity (HC2-S3, Rotronic, Switzerland) were measured at four levels (i.e., 10, 30, 50, and 70 m) of the tower, which is surrounded by typical desert grassland. The site is characterized by a middle temperate zone and semi-arid continental climate. During the study period, the daily air temperature at the 2 m height ranged between −27.3 and 33.9^∘, with an average value of 6.3^∘ (Fig. 2a). Surface-level air pressure has an inverse relation with air temperature, with an average value of 862.9 hPa (Fig. 2b). In addition, the daily average relative humidity at the 2 m height maintains a level of 41.02 % and fluctuates back and forth. The average wind speed at the 70 m height is 7.6 m s⁻¹. The daily averaged wind speed in spring and autumn occasionally exceeds the level of 10.0 m s⁻¹, indicating that the site has sufficient wind resources in these two seasons (Fig. 2c). The predominant wind direction was southwesterly and northwesterly during the whole observation period (Figs. 2d and 3).

Figure 1 (a) The observation site marked as a black spot in the Inner Mongolia Autonomous Region of China (the red line indicates the border of the Inner Mongolia Autonomous Region); (b) terrain elevation map of the 28 km × 28 km grid; (c) Google satellite historical imagery of the 28-km × 28-km grid (from © Google Maps 2021); the red dots indicate wind turbines.

Figure 2 Panels (a), (b), (c), (d), and (e) represent 2 m daily average temperature (^∘), surface-level daily average pressure (hPa), 70 m daily and monthly average wind speed (m s⁻¹), 70 m daily average wind direction (^∘), and the average 2 m daily relative humidity (%), respectively.

Figure 3 The distribution of the high-altitude wind speed and direction rose diagram of the wind measurement tower on the Ningyuanbailiutu wind farm from September 2018 to August 2020. Panels (a), (b), (c), and (d) represent 10, 30, 50, and 70 m in altitude, respectively.

2.2 Methods

2.2.1 Kernel, Weibull, and Rayleigh distributions

The kernel density estimator is the estimated PDF of a random variable. For any real values of v, the formula for the kernel density estimator is given by:

$$\begin{array}{}\text{(1)}& f_h\left(v\right)={\displaystyle \frac{\mathrm{1}}{nh}}\sum _{i=\mathrm{1}}^{n}K\left({\displaystyle \frac{v-v_i}{h}}\right),\end{array}$$

where v₁, v₂, …, v_i are random wind samples from an unknown distribution, n is the sample size, K(⋅) is the kernel smoothing function, and h is the bandwidth.

The PDF of the Weibull distribution is given by:

$$\begin{array}{}\text{(2)}& f_{\mathrm{W}}\left(v\right)=\left({\displaystyle \frac{k}{c}}\right){\left({\displaystyle \frac{v}{c}}\right)}^{k-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{v}{c}}\right)}^k\right].\end{array}$$

The Rayleigh model is a special and simplified case of the Weibull model. It is obtained when the shape factor k of the Weibull model is assumed to be equal to 2.

The maximum likelihood estimation method is a mathematical expression recognized as a likelihood function of the wind speed data in a time series format. In this method, many numerical iterations can be required to determine the k and c parameters of the Weibull function. The parameter estimation formula of the maximum likelihood method is as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }k={\left({\displaystyle \frac{\sum _{\mathrm{1}}^n{v}_i^k\ln \left(v_i\right)}{\sum _{\mathrm{1}}^n{v}_i^k}}-{\displaystyle \frac{\sum _{\mathrm{1}}^n\ln \left(v_i\right)}{n}}\right)}^{-\mathrm{1}},\text{(4)}& {\displaystyle }c={\left({\displaystyle \frac{\mathrm{1}}{n}}\sum _{\mathrm{1}}^n{v}_i^k\right)}^{\mathrm{1}/k}.\end{array}$$

The average value and standard deviation of the wind speed can be obtained from the following formulas, respectively:

$$\begin{array}{}\text{(5)}& {\displaystyle }{v}_{\mathrm{m}}={\displaystyle \frac{\mathrm{1}}{n}}\left[\sum _{i=\mathrm{1}}^n{v}_i\right],\text{(6)}& {\displaystyle }\mathit{\sigma }={\left[{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n{\left(v_i-v_{\mathrm{m}}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Alternatively, the mean wind speed can be determined from:

$$\begin{array}{}\text{(7)}& v_{\mathrm{m}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}vf\left(v\right)\mathrm{d}v,\end{array}$$

if the PDF is known.

If Eq. (7) is solved together with Eq. (2) making the substitution of $\mathit{\xi }=(v/c{)}^k$ for v, the following is obtained for the mean wind speed:

$$\begin{array}{}\text{(8)}& v_{\mathrm{m}}=c\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{k}}\right).\end{array}$$

Note that the gamma function has the properties of $\mathrm{\Gamma }\left(x\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\xi }}^{x-\mathrm{1}}\exp (-\mathit{\xi })\mathrm{d}\mathit{\xi }$ and $\mathrm{\Gamma }(\mathrm{1}+x)=x\mathrm{\Gamma }\left(x\right)$.

2.2.2 Power density

The mean power density for the kernel smoothing function becomes:

$$\begin{array}{}\text{(9)}& P_{\mathrm{K}}=\sum _{i=\mathrm{1}}^n\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{v}_{\mathrm{m},i}^{\mathrm{3}}{f}_h\left(v_i\right)\right].\end{array}$$

The mean power density for the Weibull function becomes:

$$\begin{array}{}\text{(10)}& P_{\mathrm{W}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{c}^{\mathrm{3}}\mathrm{\Gamma }\left[\mathrm{1}+{\displaystyle \frac{\mathrm{3}}{k}}\right].\end{array}$$

The mean power density for the Rayleigh model is found to be:

$$\begin{array}{}\text{(11)}& P_{\mathrm{R}}={\displaystyle \frac{\mathrm{3}}{\mathit{\pi }}}\mathit{\rho }{v}_{\mathrm{m}}^{\mathrm{3}},\end{array}$$

where ρ is the air density.

2.2.3 Weibull parameters

The relationship between scale factor c and height can be expressed as follows:

$$\begin{array}{}\text{(12)}& c/c_{\mathrm{10}}=(z/\mathrm{10}{)}^{\mathit{\alpha }}.\end{array}$$

Here c₁₀ represents the scale factor at 10 m height, z represents the height, and α represents the power exponent parameter to be estimated.

The relationship between scale factor k and height can be expressed as follows:

$$\begin{array}{}\text{(13)}& k=a(z/\mathrm{10}{)}^{\mathrm{2}}+b^{\prime }(z/\mathrm{10})+d,\end{array}$$

where a, b^′, and d are unknown parameters to be fitted to the quadratic function.

In addition, as shown in Fig. 6c, Justus and Mikhail (1976) gave the following formula for the shape factor k with height:

$$\begin{array}{}\text{(14)}& k_{\mathrm{10}}/k=\mathrm{1}+b_{\mathrm{10}}\ln (z/\mathrm{10}),\end{array}$$

where k₁₀ is the shape factor at a reference height of 10 m. At a reference height of 10 m, b=b₁₀ is just some constant, whose value can be determined by a least squares fit of relation (14) to the data.

2.2.4 Surface roughness length

When the wind speed at three or more heights is measured, the roughness length calculated by the least square regression (Archer and Jacobson, 2003; Archer, 2005; Golbazi and Archer, 2019) is:

$$\begin{array}{}\text{(15)}& \ln \left(z_{\mathrm{0}}\right)={\displaystyle \frac{\begin{array}{c}U\left(z_{\mathrm{R}}\right)\left\{\sum {\left[\ln \left(z_i\right)\right]}^{\mathrm{2}}-\ln \left(z_{\mathrm{R}}\right)\sum \ln \left(z_i\right)\right\}\\ -\ln \left(z_{\mathrm{R}}\right)\sum \left[U_i\ln \left(\frac{z_i}{z_{\mathrm{R}}}\right)\right]\end{array}}{\begin{array}{c}U\left(z_{\mathrm{R}}\right)\sum \ln \left(z_i\right)-\sum \left[U_i\ln \left(\frac{z_i}{z_{\mathrm{R}}}\right)\right]\\ -N\cdot U\left(z_{\mathrm{R}}\right)\ln \left(z_{\mathrm{R}}\right)\end{array}}},\end{array}$$

where z_R is the reference height, z_i is the height of the other three layers, N=4 representing four vertical layers, and U_i is the wind speed corresponding to the height of the other three layers. In most cases, it is a purely mathematical statistical method, and therefore this simple mathematical method does not require a physical explanation for roughness estimation.

In addition, the aforementioned method is obtained from the logarithmic wind speed profile, which is a typical form of wind speed profile under neutral stratification. A calculation of the wind speed at other altitudes under the reference altitude can be obtained from the following formula (Golbazi and Archer, 2019; Archer and Jacobson, 2003):

$$\begin{array}{}\text{(16)}& U\left(z\right)=U\left(z_{\mathrm{R}}\right){\displaystyle \frac{\log \left(\frac{z}{z_{\mathrm{0}}}\right)}{\log \left(\frac{z_{\mathrm{R}}}{z_{\mathrm{0}}}\right)}},\end{array}$$

where z₀ is the estimated surface roughness length, assuming that the friction speed near the ground does not change with height.

3 Results

3.1 Comparisons of kernel, Weibull, and Rayleigh models

The monthly, seasonal, and annual average wind speed values and standard deviations calculated using Eqs. (5) and (6) for the available time series data are shown in Table 1. It can be seen from Table 1 that the highest average wind speeds occurred in May and December 2019 and in May 2020, and the lowest average wind speeds occurred in February and August 2019. Over about 2 years, it was found that the average wind speed in the spring of 2019 and 2020 was higher, and the average wind speed in the summer of 2019 and 2020 was lower. During the entire study period, the average wind speed values at 10, 30, 50, and 70 m were 6.0, 6.8, 7.2, and 7.6 m s⁻¹, respectively, which also shows that the wind speed value increases with an increase in altitude.

Table 1 Calculated monthly, seasonal, and annual distribution parameters based on the time series wind speed data measured every 15 min from Damaoqi Wind Tower. (The red shading and blue shading represent the larger and smaller values in the table, respectively. The darker the color, the more extreme the value.)

Figure 4 Frequency density histogram of wind speed at 70 m height from autumn 2018 to summer 2020; the probability density curve obtained by fitting the Weibull, kernel, and Rayleigh distribution functions to the original data. Panels (a), (b), (c), (d), (e), (f), and (h) represent each season, respectively.

Figure 4 shows the frequency density histogram of the wind speed at 70 m for about 2 years and the probability density curves of the Weibull, kernel, and Rayleigh distributions. First of all, it is obvious from the frequency histogram that the wind speed at 70 m fluctuated drastically in the autumn of 2018, spring of 2019, and summer of 2020. This conclusion can also be confirmed from the data in Table 1. The shape factor k values of these three specific seasons are 2.18, 2.13, and 2.11, respectively, which are slightly higher than the shape factor k value of the Rayleigh distribution. In combination with Table 1, it is also found that the higher the value of the scale factor c, the smoother the three specific probability distribution curves. By contrast, as shown in Fig. 4d, its three specific probability density curves are very sharp.

Although the kernel distribution also has specific parameters to control its probability density curve, it does not have the general form of wind speed distribution. Moreover, the k value of the Weibull distribution is ∼ 2. To select the specific wind speed distribution form suitable for the Ningyuanbailiutu site, therefore, the model prediction accuracies of the Weibull distribution and the Rayleigh distribution for average wind power need to be compared further.

Figure 5 (a) The average power density values calculated by the four distributions: the frequency distribution of the original wind speed time series data, the kernel distribution, the Weibull distribution, and the Rayleigh distribution, for each month, each season, and the total period of about 2 years. The dotted red line represents the actual grid-connected average power density when the fan blade length is assumed to be 41.5 m. Panels (b), (c), and (d) represent the residual and 95 % confidence interval under the three specific probability model distributions, respectively. The unit of the legend is W m⁻². The total average power density on the right y axis in (a) is the average power density of the grid-connected wind turbine with a radius of 41.5 m represented by the red dots in Fig. 1c. The blue left y axis in panels (b), (c), and (d) matches the shaded part formed by the blue line, which represents the average power density; the red right y axis matches the red rectangle, which represents the residual error between the model value and the actual value.

The mean power densities calculated from the measured probability density distributions and those obtained from the models are shown in Fig. 5. The mean power density shows significant monthly and seasonal variation. The minimum average power density appeared in August 2019 and was only 214.9 W m⁻². In addition, smaller mean wind power densities appeared in July and September 2019 and January, July, and August 2020, which were generally lower than 350.0 W m⁻². Generally, the maximum value of monthly mean wind power density reached 862.4 W m⁻² in May 2019, and the seasonal mean wind power density peaked in spring 2020.

The differences between the kernel distribution, Weibull distribution, and Rayleigh distribution are explored when calculating the average wind power density and the frequency distribution using the original wind speed data. The 2-year mean absolute percentage error (MAPE) values in calculating the mean power density using the kernel, Weibull, and Rayleigh functions are 1.17 %, 1.05 %, and 4.20 %, respectively. The RMSEs of the kernel distribution, Weibull distribution, and Rayleigh distribution are 45.8, 60.5, and 875.3 W m⁻², respectively. The Weibull and kernel models return smaller error values in calculating the mean power density compared to the Rayleigh model. The mean power density is estimated by the Rayleigh model to have a very large absolute error value of 83.1 W m⁻² in December 2019. On the other hand, the highest absolute error value occurs in May 2019 with 21.3 W m⁻² for the Weibull model.

Analysis of residual error and average percentage error suggests that the average wind power density estimated by the Weibull distribution with specific parameter control is very similar to the kernel distribution, which is closest to the original wind frequency distribution (Fig. 5c). The lower limit of the 95 % prediction interval is each predicted value minus 1.96 standard deviations, and the upper limit is each predicted value plus 1.96 standard deviations (Fig. 5b–d). This suggests that the interval applicability of the three specific distribution models is good.

3.2 Vertical characteristics of Weibull parameters

Figure 6 shows the characteristic variation of the scale factor c and the shape factor k with height estimated from the Weibull distribution for original wind speed data during the study period, exhibiting power exponential and quadratic function variations, respectively.

Figure 6 (a) The characteristic variation of scale factor c with height based on Eq. (12). The characteristic variation of shape factor k with height based on (b) Eq. (13) and (c) Eq. (14).

Table 2 The values of various parameters in different time periods and the corresponding RMSEs under least squares formula fitting. (The red shading and blue shading represent the larger and smaller values in the table, respectively. The darker the color, the more extreme the value.)

Table 2 gives in detail the values of α, a, b^′, d, and b₁₀ obtained by the least squares fitting method for each month, each season, and all time periods, as well as the corresponding RMSEs obtained from the formula. When using the power exponent formula (12) to fit the relationship between the scale factor c and the height, the RMSE has the smallest values in January 2019 and July and August 2020. However, in December 2019, January 2020, and February 2020 it has the largest values. This shows that formula (12) has a better fitting effect in the winter of 2018, and a poor fitting effect in the winter of 2019. Justus and Mikhail (1976) found that the mean value of α was 0.23. In the present study, the mean value of a for each month over the whole 2 years is 0.117, and the corresponding standard deviation is 0.016.

Figure 6b and c indicate that the two different formula forms have a good fitting relationship for shape factor k and height. The RMSEs of Table 2 also suggest that the effect of the quadratic function fitting is better than the logarithmic reciprocal function of Justus and Mikhail (1976). The RMSE of the quadratic function fitted to all data for 2 whole years is 0.0078, but the RMSE of the logarithmic reciprocal function is 0.0214, which is close to a multiple of 1 : 3. Both of these two types of formula are basically applicable only to heights below 100 m. In addition, from a comparison of Fig. 6b and c, it can be seen that there will be different trends in the change in the k value with height, and the increasing or decreasing speed of the k value in the form of a quadratic function will be higher than that found by Justus and Mikhail (1976) when the height is greater than 70 m. This different trend will lead to large errors in estimating wind energy resources above 70 m.

3.3 Spatial–temporal variations in surface roughness length

The shape of the wind profile is greatly affected by the surface roughness in the direction of the incoming flow. Thus, surface roughness is a key element in wind energy resource evaluation and forecasting models. In calculating aerodynamic roughness, especially in practical applications, the least squares approximation of the logarithmic profile equation to the measured wind speed profile method has been widely used, referred to as the “logarithmic profile method”.

Figure 7 The average and median values of estimated roughness in each month for the total period.

After calculating the 15 min continuous wind speed data using the above method, quality control of the data is carried out. In this study, we eliminated wind speeds greater than or equal to 6.0 m s⁻¹ at 50 m, and the estimated abnormal roughness data are infinitely large or infinitely small. Figure 7 shows that both the average and median monthly roughness lengths in January, February, and March 2019 are significantly less than those in August, September, and October 2019. The largest value of median roughness was close to 0.19 m in October 2019, and the maximum value of average roughness was approximately equal to 0.27 m. In June 2020, the median and average roughness values reached 0.18 and 0.25 m, respectively. The minimum value of median roughness was about 0.10 m in January 2019, and the smallest value of average roughness value was about 0.20 m in January 2020.

In addition, the median and average roughness length were lowest at about 0.12 and 0.22 m in the winter of 2018 and 2019, while the highest were about 0.15 and 0.25 m in the autumn of 2019. It is notable that the roughness length steadily increases from winter to autumn. In short, this suggests that the grassland vegetation in the site area has an obvious wilting period, growing season, and lush period. Compared with the average roughness length, the representative roughness length of the area fitted the median value more closely.

According to the Davenport land type roughness classification (Davenport et al., 2000) and the summary of roughness length over the wind tower sites and the corresponding land types (Li et al., 2021), in the case of land types with less vegetation and cropland, the roughness length is generally estimated to be a slightly rough open area of about 0.10 m. The area we studied belongs to the grassland vegetation type, and the roughness estimate should be around 0.13 m, and will not be classified as rough; that is, the roughness length is as high as 0.25 m. In addition, in a study (Golbazi and Archer, 2019) on the estimation of sea surface roughness length in coastal waters, it is mentioned that the statistical method uses a single constant value of z₀ in the representative area, and the median value can be worth recommending.

Figure 8 The average and median values of estimated roughness in 12 different directions of incoming flow.

Figure 8 shows the variation in the estimated roughness length in 12 different incoming wind directions. When the wind direction is 120 or 240^∘, the estimated roughness length is highest, and the median value and average value are about 0.23 and 0.30 m, respectively. Secondly, when the wind direction is 30 or 300^∘, the estimated roughness length is lowest, and the median value and average value are about 0.08 and 0.18 m, respectively. Therefore, between the highest and lowest estimated roughness lengths, there is a specific trend of increasing or decreasing. This phenomenon can be explained in conjunction with Figs. 1 and 3. There is a hillside to the west of the wind tower. Therefore, when the incoming wind direction is 120 or 240^∘, it is on the windward side or leeward side, respectively, of the wind-measuring tower. In this way, there will be a pressure difference, which will increase friction loss and increase the estimate of the effective roughness length. When the incoming wind direction is 30 or 330^∘, it is found that the wind passing through the wind measurement tower will not be greatly affected by the terrain. The terrain is relatively flat, and the estimated roughness length is close to the normal value of 0.10 m. In addition, in the plot of roughness length estimation with wind direction, there are obviously more data points in the wind directions from 180 to 330^∘ than in the other wind directions. The 240^∘ wind direction has the most data points, which also shows that the site has a southwesterly wind blowing all year round.

3.4 Extrapolation of the average wind power density

With the scale factor c changing with height in the form of a power function, and shape factor k changing with height in the form of a quadratic function, the scale factor c and the shape factor k at 75, 80, 85, 90, 95, and 100 m are calculated. Then the average wind power density (Fig. 9b) is calculated for each month, each season, and the whole time period from formula (10). On the other hand, when studying the roughness length parameter in the previous section, we assume that the roughness length calculated from the four-layer height is dynamic. Then through the logarithmic form of formula (16), we can calculate the wind speed values at 75, 80, 85, 90, 95, and 100 m every 15 min. Finally, the “reference average power density” (Fig. 9a) at six specific heights can also be obtained.

Figure 9 The logarithmic model (a), power-law model (b), and their residuals (c) estimate the average wind energy density at six specific heights (75, 80, 85, 90, 95, and 100 m) of the wind measurement tower in each month, for every season, and over the whole period.

Both the power-law model and the logarithmic model can estimate the average wind power density of six specific heights, and as seen in Fig. 9 the values estimated by the two methods show greater differences in autumn and winter, and smaller differences in spring and summer. In addition, the two different models both show that the average power density values are largest in spring and smallest in summer. Although the average power density values increased with height over the whole experiment, the gradient of the increase in average power density values with height is largest in autumn and winter, and smallest in spring and summer. Figure 9c shows that relative to the power-law model, the average power density of the logarithmic model extrapolated at 70–100 m is smaller in the winter of 2018 and in July and August of 2020, while it is larger in other study periods. Generally speaking, the difference between the estimated average power density values is very small. However, the data and methods used in the estimation of the two models are different. The result of this estimation gives us important guidance for studying two Weibull parameters, namely, the scale factor c and the shape factor k, and the surface roughness length parameter.

4 Discussion

4.1 Applicability of Weibull and Rayleigh models

There are various statistical distribution functions for describing and analyzing wind data, including normal, lognormal, Rayleigh, and Weibull probability distributions (Fagbenle et al., 2011; Ozerdem and Turkeli, 2003). It has been found that the Weibull and Rayleigh distributions are the most accurate and adequate in wind analysis as well as in interpreting the actual wind speed data and in predicting the characteristics of the prevailing wind profile. A kernel distribution is a nonparametric representation of the PDF of a random variable. A kernel distribution is defined by a smoothing function and a bandwidth value, which control the smoothness of the resulting density curve (Kafadar et al., 1999). In fact, some scholars have used the probability density distribution of wind speed to compare the advantages and disadvantages of the Weibull distribution and the Rayleigh distribution (Celik, 2004; Pishgar-Komleh et al., 2015). In our present work, the kernel function exhibits the feature of the smooth function, and is also closer to the actual frequency distribution (Fig. 4), which can be used to fit the original wind speed data. Therefore, the kernel function is employed as a medium for comparing the pros and cons of the Weibull and Rayleigh functions in the desert steppe area.

Celik (2003) employed the Weibull function to analyze wind power density in six different regions around the world, and the average percentage errors obtained were relatively low. The reason may be that the scale factor c representing the average wind and the shape factor k are relatively small, the original wind speed is relatively stable and fluctuates little, and also the Weibull function has strong applicability. Celik (2004, 2011) used the Weibull function and Rayleigh function to calculate wind power density in Turkey, showing the average annual relative errors of < 8 % and 37 %, respectively. Pishgar-Komleh et al. (2015), on the other hand, reported annual average errors of 55.00 % for both the Weibull function and the Rayleigh function. These values vary largely from our present results, due to the different applicability of the two specific distribution functions in different regions. Therefore, when developing new wind farms, it is extremely important and necessary to compare the applicability of various wind speed distribution functions in the local area. In general, we found that the Weibull distribution is applicable for depicting the wind speed distribution at the Ningyuanbailiutu site in northern China.

In addition, although the mean wind power density calculated in this study is in good agreement with the actual grid-connected average power density (Fig. 5a), there is a significant difference in these two values. This is because the wind turbines are not always connected to the grid, due to failures, or to other wind turbines outside the range of the wind measurement tower. As a result, the wind measured by a single wind tower will underestimate the wind speed of other wind turbines. The limitation of this study is that data from only a single wind tower were derived. For future wind energy density estimations, it is worth collecting data from more wind towers to obtain a more realistic wind resource distribution in the study area.

4.2 The complexity of Weibull parameters and surface roughness

Figures 6, 7, 8, and Table 2 show the spatiotemporal variability of scale factor c, shape factor k, and surface roughness z₀, which can be attributed to the following three aspects: (a) the type of surface land and meteorological conditions (Golbazi and Archer, 2019); (b) the uncertainty of Weibull parameters calculated using maximum likelihood method (Mohammadi et al., 2016) and uncertainty of the roughness length calculated using statistical mathematical methods (Kim and Lim, 2017); (c) the limitations of extrapolating high-level Weibull parameters (Justus and Mikhail, 1976). The shape factor k varies with height and exhibits not only the form of the reciprocal of the logarithmic function but also the form of the quadratic function. From the RMSE in Table 2, it can be seen that the quadratic function is the most suitable for this study area. Therefore, to use the Weibull function to evaluate the high-level wind speed distribution in a specific research area in the future, we should consider high-density observations with more fitting methods so as to obtain the best functional form of the Weibull parameter varying with height. Table 3 reviews scale factor c, shape factor k, surface roughness z₀, and yearly mean absolute percentage error over different topographies, showing obvious regional differences, due to various climate and topography contexts with different methods. For example, Pishgar-Komleh et al. (2015) used a constant surface roughness value of 0.14 to extrapolate wind speed, ignoring the dynamic changes in surface roughness throughout the year. The calculated annual mean absolute percentage error is much higher than the value calculated by the dynamic surface roughness in this study. Therefore, we should take dynamic roughness into account based on a reliable and accurate topographic map, rather than assuming surface roughness as a constant.

Table 3 Review of scale factor c, shape factor k, surface roughness z₀, and yearly mean absolute percentage error (MAPE) over different topographies.

The differences and uncertainties between the logarithmic and power-law models can also be seen in Fig. 9. Its uncertainty is manifested in the absence of verification of actual high-level wind speeds. The reason for the difference is that the shape factor k in the winter of 2018 and in July and August of 2020 shows a decreasing trend with height, and the average wind energy density is inversely proportional to the shape factor k, according to formula (10). Therefore, the limitation of this paper is that the extrapolated results need to be further confirmed by future encrypted observations of high-level wind speed data. In addition, the encouraging news is that this gives us two reliable options for future extrapolation of high-level wind energy density: In the case of non-Weibull winds, only the logarithmic model can be considered, and in the case of Weibull winds, both the logarithmic model and the power-law model are good choices.

5 Conclusions

The present work investigated the scale factor c and the shape factor k that affect the Weibull distribution of wind speed, by directly estimating the energy potential of the wind speed resource at four different heights, and the surface roughness length parameter that directly affects the shape and law of the wind profile. The main conclusions are:

The 2-year mean absolute percentage error values in calculating the mean power density using the kernel, Weibull, and Rayleigh functions are 1.17 %, 1.05 %, and 4.20 %, respectively. The Weibull wind speed distribution model is the most suitable wind speed distribution model for the Ningyuanbailiutu site. The scale factor c increases with an increase in height, showing an obvious form of power function. The shape factor k increases or decreases with height and has two different forms, which are the reciprocal of the quadratic function and the logarithmic function. For further determination of the changes in form factor with height, it will be necessary in future to set up intensive observations for heights above 70 m and below 100 m.

When estimating the surface roughness length, the median value is selected as the representative value of the surface roughness length. This is based not only on recognition of actual previous research, but also on confirmation of actual grassland vegetation types. Although the statistically calculated z₀ does not have a proper physical explanation, it gives the most accurate wind speed estimate at the required height. The estimated roughness length varies with the seasons of the grassland vegetation at the site. The estimated roughness lengths of the wilting period, growing season, and lush period are about 0.12, 0.13, and 0.15 m, respectively. The estimated surface roughness length will be affected by the windward and leeward sides. When the wind flows across the hillside, there will be a pressure difference, which will increase the friction loss and increase the estimated effective roughness length. The prevailing wind direction at this site is 240^∘, which happens to be the direction of the windward side of the site. The estimated roughness length is about 0.23 m. Finally, the power-law model and the logarithmic model were employed to estimate the average power density values at 75, 80, 8, 90, 95, and 100 m. The two models show greater differences in autumn and winter, and smaller differences in spring and summer. The gradient of the increase in average power density values with height is largest in autumn and winter, and smallest in spring and summer.

In general, against a carbon-neutral background, the determination of the potential for economical and clean wind energy resources is an important scientific issue in the development of renewable energy worldwide. Our research has determined the possible relationship between Weibull natural wind mesoscale parameter c and shape factor k with height under the conditions of a desert steppe terrain in northern China; this has great potential in wind power generation, but there is a lack of comprehensive investigations into key parameters for estimating wind power density from tower data. In the present study, we gained an enhanced understanding of the seasonal changes in the surface roughness of the desert grassland and the changes in the incoming wind direction. Our findings also have important implications for the assessment of wind energy resources for the establishment of new wind farms in areas experiencing varied desert steppe terrains throughout the world.