1 Introduction

To reduce greenhouse gas emissions, harnessing renewable and sustainable energy resources, such as wave energy, wind energy, and solar energy, has garnered increasing interest (Wu et al., 2021; Yang et al., 2021; Cheng et al., 2022a, 2022b; Abu-Rayash and Dincer, 2023). Offshore wind energy is a widely recognized renewable energy source because of its abundant resource potential, mature technology, and relatively low cost. Over the past decade, the field has witnessed considerable attention and development, with its expansion being anticipated to rise in the coming years (Costoya et al., 2022; Martinez and Iglesias, 2022; Zhuo et al., 2022). Bilgili and Alphan (2022) reported that the total capacity of offshore wind power will increase to approximately 382 GW by 2030 and 2002 GW by 2050. Note that an in-depth wind energy appraisal at the intended hub height of wind turbines is imperative for selecting a wind turbine model and contributes to the economic feasibility analysis of offshore wind power projects (Jung and Schindler, 2017; Nie and Li, 2018; Wang et al., 2022a, 2022b; Wu et al., 2022a).

Over the last few decades, different resource assessment factors, such as wind power density (WPD), the occurrence rate of exploitable wind energy, and the variability index, have been widely used in the wind energy industry to measure wind resource potential (Zheng et al., 2013; Langodan et al., 2016; Allouhi et al., 2017; Wen et al., 2021). Among them, WPD is the most crucial metric because most other factors can be derived from this wind energy variable. Due to the cubic relationship between wind speed and WPD, high-quality long-term wind speed data at corresponding hub heights are typically required to compute the energy production (EP), capacity factor (CF), and other evaluation parameters related to wind turbine production (Nie and Li, 2018; Yang et al., 2023). Wind speed data can be derived from various sources such as meteorological stations or buoys (Jiang et al., 2013; Liang et al., 2020; Dayal et al., 2021), publicly accessible reanalysis wind field datasets (Zheng et al., 2013, 2018; Wen et al., 2021), and numerical models (Langodan et al., 2016; Dong et al., 2019; Li et al., 2022). However, they offer wind data at fixed hub heights, which warrants the extrapolation of the original wind speeds to intended hub heights. The most commonly utilized method for extrapolation is the Hellman power law profile method (Rehman and Al-Abbadi, 2005; Jung and Schindler, 2017; Gualtieri, 2019; Kavari et al., 2019; Han et al., 2022; Yan et al., 2022), in which the wind shear coefficient is the sole parameter of this method. The wind shear coefficient represents the speed gradient in the vertical direction and determines the quality of the extrapolated wind speed (Ray et al., 2006). As a consequence, the joint behavior of wind speed and wind shear coefficient (WS-SC) is crucial for height-dependent wind resource characteristics analysis.

So far, many parametric models, such as Weibull, lognormal, Rayleigh, and Gamma, have been used for constructing the probability distribution function (PDF) of wind speed, as demonstrated by Carta et al. (2009) and Jung and Schindler (2019b). However, because of the combination of multiple weather systems, the probability distribution of wind speed may exhibit multimodal structures. Traditional single-mode distributions may thus be insufficient in providing a satisfactory description of the complex patterns in the actual wind field. The finite mixture model offers a practical solution in which a linear combination of several aforementioned parametric models is used to establish the distribution. In the past decades, several mixture models, including two-component mixture Weibull distribution, Weibull singly truncated normal mixture distribution, and Weibull lognormal mixture distribution, have been set up and employed to investigate wind speed data (Carta and Ramírez, 2007; Akpinar and Akpinar, 2009; Chang, 2011; Kollu et al., 2012; dos Santos et al., 2021).

For the other important parameter in height-dependent wind resource assessment, the wind shear coefficient is usually assumed to be a constant of 0.143 for an onshore wind energy analysis and 0.11 for an offshore wind power investigation (Hsu et al., 1994; Shu et al., 2015; Kruyt et al., 2017; Mahmoodi et al., 2020; Wen et al., 2022). However, the value of this parameter can be affected by wind speed, local surface roughness, temperature, pressure, atmospheric stability, season, and terrain type (Irwin, 1979). The magnitude of this variable also depends on the nonlinear wind-wave interaction and the heat capacity of water for offshore sites (Tambke et al., 2005). Therefore, the assumption of constants does not hold in general. The effect of the wind shear coefficient on wind resource analysis has been examined and has recently gained eminence (Rehman and Al-Abbadi, 2005, 2007; Fırtın et al., 2011; Gualtieri, 2016; Werapun et al., 2017; Li et al., 2018; Emeksiz and Cetin, 2019; Kim et al., 2021). Werapun et al. (2017) observed that ignoring the variation in the wind shear coefficient can result in errors of up to 35% in EP estimation at a specific hub height, while Fırtın et al. (2011) showed that this error can be as high as 49.6%. Yan et al. (2022) examined offshore Lidar measurements and disclosed that wind shear coefficient values are subject to seasonal variations. Interestingly, the probability distribution of the wind shear coefficient also shows a multimodal structure. In this regard, a mixture model is crucial for statistical modeling (Corrigan and Matthiesen, 2017). The results of these studies underscore the significance of accurately accounting for wind shear coefficient variation in wind resource characteristic investigation to acquire reliable estimates of the energy potential.

To accurately evaluate wind resources, a suitable bivariate distribution model is required to fully capture the joint statistical properties of the WS-SC. The copula-based method has emerged as an effective approach for building bivariate distribution models over the past decades (Dong et al., 2015; Lucas and Guedes Soares, 2015; Huang and Dong, 2021b; Song et al., 2022). This method involves fitting marginal distributions separately using univariate distribution models and describing the dependency of variables using a copula function. Han et al. (2022) applied the Dagum-normal distribution and kernel density estimation for the WS-SC and linked them via a nonparametric copula function. They suggested that this bivariate distribution model can boost the accuracy of site selection and further enhance the economic benefits of wind power projects. By analyzing the empirical probability density of WS-SC, Jung and Schindler (2017) observed that the dependency between these variables changes at 9.0 m s⁻¹. As a consequence, they separated the original bivariate sample into two parts. For each subsample, the four-parameter Johnson SB and Dagum distributions were utilized as the marginal distribution models for WS-SC, respectively, and the corresponding dependence structure was described by a Gaussian copula function. However, it should be noted that although the multimodal characteristics of bivariate data are considered, using only a constant of wind speed to separate the modes involved in the data is not logical. Therefore, the constructed models cannot exhibit satisfactory performances at all sites.

Recently, the mixture copula method has become known as a promising statistical modeling approach to characterizing bivariate data showing multimodal features. This method has found increasing application in establishing multivariate joint distribution models for wave spectrum parameters in combined sea states. Lin et al. (2020) categorized the wave data into four components based on the values of wind sea and swell energy and then utilized a combination of copula functions to model the statistical patterns of wave height and period. Huang and Dong (2021a), in contrast, put forward a method for constructing a bivariate distribution model using total data. The number of mixture components was determined according to the Bayesian information criterion (BIC). Both methods have shown superior performance in describing the statistical patterns of ocean environment parameters compared with the commonly used single copula model and conditional modeling approach. More recently, Yang et al. (2023) implemented a mixture bivariate distribution model to describe the probabilistic structures of wind speed and air density, realizing enhanced accuracy in wind energy evaluation.

This work endeavors to develop a joint statistical model of WS-SC using the mixture copula method. The objective of this study is to facilitate offshore wind power development by conducting an in-depth height-dependent wind resource characteristic analysis. The major novelty and scientific contributions of this work are as follows:

1) A mixture bivariate parametric model was developed to model the joint distribution of WS-SC under mixed atmospheric conditions. The model parameters were estimated using the Expectation Maximization (EM) algorithm, and BIC was adopted to determine the optimal number of components.

2) The multimodal structures shown in the joint probability distribution of WS-SC under mixed atmospheric conditions can be captured well by the proposed model.

3) The wind resource characteristics at arbitrary heights can be computed after one parametrization of the proposed mixture model. Based on this, the temporal and spatial distribution characteristics of height-dependent wind resources were investigated by considering the variation in the wind shear coefficient.

4) The influence of wind shear variation on wind resource characteristics was quantified and analyzed.

The rest of this paper is structured as follows. Section 2 presents the methods for building the mixture bivariate statistical model, performing height-dependent wind energy evaluation, and estimating the effect of wind shear coefficient variation. Section 3 presents the data used in this work and their analysis. In Section 4 the proposed mixture copula model is used to build the joint statistical model of WS-SC. The height-dependent wind resource characteristics are analyzed to provide insights into the potential for wind energy generation at various heights. This is followed by a quantitative analysis of the effect of wind shear variation on wind resource characteristics. Section 5 summarizes the major findings and their implications for offshore wind power development.

2 Methodology 2.1 Joint Statistical Model for Wind Speed and Wind Shear 2.1.1 Copula theory

The copula theory proposed by Sklar (1959) has emerged as a potent method for constructing multivariate distribution models, allowing the joint statistical modeling of several dependent random variables by integration of multiple marginal distributions using a copula function, which is a measure of dependence. Here, we adopted X and Y as wind speed (V) and wind shear coefficient (E), respectively. The bivariate distribution of the WS-SC can be expressed as

$ {F_{XY}}(x, y) = {C_{XY}}\left[ {{F_X}(x), {F_Y}(y)} \right], $

(1)

where F_XY(x, y) is the bivariate cumulative distribution function (CDF) of WS-SC, F_X(x) and F_Y(y) are the marginal CDFs of WS-SC, respectively, and C_XY[·] is a copula function. The corresponding PDF is given by

$ {f_{XY}}(x, y) = {c_{XY}}\left[ {{F_X}(x), {F_Y}(y)} \right]{f_X}(x){f_Y}(y), $

(2)

where f_X(x) and f_Y(y) are the marginal PDFs of the WS-SC, respectively, and c_XY[·] is the copula density function

$ {c_{XY}}\left[ {{F_X}(x), {F_Y}(y)} \right] = \frac{{{\partial ^2}{C_{XY}}\left[ {{F_X}(x), {F_Y}(y)} \right]}}{{\partial {F_X}(x)\partial {F_Y}(y)}}. $

(3)

Different families of copula functions have been developed for statistical modeling (Guedes Soares and Carvalho, 2012; Dong et al., 2015; Vanem, 2016; Bai et al., 2020; Lin et al., 2020; Zhai et al., 2021; Gong et al., 2022). Among them, the Gaussian copula function is preferred (Lin et al., 2020; Huang and Dong, 2021a), which has the advantage of flexibility in capturing complex dependencies while keeping computational efficiency. Thus, in this work, we used the Gaussian copula function to jointly model wind speed and wind shear coefficient

$\begin{array}{l} {F_{XY}}(x, y) = \varPhi \left\{ {{\varPhi ^{ - 1}}\left[ {{F_X}(x)} \right], {\varPhi ^{ - 1}}\left[ {{F_Y}(y)} \right]} \right\} =\\ \int_{ - \infty }^{{\varPhi ^{ - 1}}\left[ {{F_X}(x)} \right]} {\int_{ - \infty }^{{\varPhi ^{ - 1}}\left[ {{F_Y}(y)} \right]} {\frac{1}{{2{\text{π }}\sqrt {1 - {\gamma ^2}} }}} } \exp \left[ { - \frac{{{u^2} - 2\gamma uv + {v^2}}}{{2(1 - {\gamma ^2})}}} \right]{\text{d}}u{\text{d}}v. \end{array} $

(4)

where Φ{·} is the CDF of the standard normal distribution and γ is the parameter of the Gaussian copula function.

2.1.2 Univariate distribution models

Statistical modeling of wind speed has been paid considerable attention, with several parametric distribution models vying for prevalence. Among them, the Weibull distribution has become a prominent contender, as shown by extensive research in the field (Carta et al., 2009). Notably, Li et al. (2021) suggested the inclusion of the Weibull distribution in a mixture model to accurately describe the distribution characteristics of wind speed. In this context, our work adopted the Weibull distribution model to describe the distribution characteristics of wind speed. The PDF and CDF of the Weibull distribution are

$ f(x|\alpha, \beta) = \frac{\alpha }{\beta }{\left({\frac{x}{\beta }} \right)^{\alpha - 1}}\exp \left[ { - {{\left({\frac{x}{\beta }} \right)}^\alpha }} \right], $

(5)

$ F(x|\alpha, \beta) = 1 - \exp \left[ { - {{\left({\frac{x}{\beta }} \right)}^\alpha }} \right], $

(6)

where α > 0 denotes the shape parameters and β > 0 is the scale parameter.

Considering the documented prevalence of the negative wind shear coefficient (Fırtın et al., 2011; Yan et al., 2022), the normal distribution was used to represent wind shear characteristics. The PDF and CDF of the normal distribution are

$ f(y|\mu, \sigma) = \frac{1}{{\sqrt {2{\text{π }}} \sigma }}\exp \left[ { - \frac{{{{(y - \mu)}^2}}}{{2{\sigma ^2}}}} \right], $

(7)

$ F(y|\mu, \sigma) = \frac{1}{{\sqrt {2{\text{π }}} \sigma }}\int_{ - \infty }^y {\exp \left[ { - \frac{{{{(y - \mu)}^2}}}{{2{\sigma ^2}}}} \right]{\text{d}}y}, $

(8)

where μ is the mean of y and σ is the standard deviation of y.

2.1.3 Mixture copula model

It has been established that WS-SC typically exhibits a multimodal univariate distribution (Corrigan and Matthiesen, 2017). As a result, the joint distribution of these two dependent random variables may exhibit a similar multimodal structure that cannot be deftly represented by a unimodal bivariate distribution model. To overcome this problem, we put forward a mixture copula distribution model that incorporates N_c components, of which the CDF is

$\begin{array}{l} {F_{XY}}(x, y|\psi) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{F_{X{Y_k}}}(x, y|{\xi _k})}\\ \qquad= \sum\limits_{k = 1}^{{N_c}} {{w_k}{C_{X{Y_k}}}\left\{ {{\varPhi ^{ - 1}}\left[ {{F_{{X_k}}}(x)} \right], {\varPhi ^{ - 1}}\left[ {{F_{{Y_k}}}(y)} \right]} \right\}}, \end{array} $

(9)

where $ \psi = {\text{\{ }}{w_1}, {\text{ }}{w_2}, \cdots, {\text{ }}{w_{{N_c}}}{\text{, }}{\xi _1}, {\text{ }}{\xi _2}, \cdots, {\text{ }}{\xi _{{N_c}}}{\text{\} }} $ denotes the model parameters, w_k refers to the weight parameter satisfying $ \sum\limits_{k = 1}^{{N_c}} {{w_k} = 1} $, ξ_k denotes the parameters of $ {F_{X{Y_k}}} $, $ {F_{X{Y_k}}} $ means the joint CDF of the kth distribution model, $ {C_{X{Y_k}}} $ refers to the copula function for the kth component, and $ {F_{{X_k}}} $ and $ {F_{{Y_k}}} $ refer to the marginal CDFs of X and Y for the kth component, respectively. The corresponding PDF f_XY

$\begin{array}{l} {f_{XY}}(x, y|\psi) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{f_{X{Y_k}}}(x, y|{\xi _k})}\\ \qquad = \sum\limits_{k = 1}^{{N_c}} {{w_k}{c_{X{Y_k}}}\left[ {{F_{{X_k}}}(x), {F_{{Y_k}}}(y)} \right]} {f_{{X_k}}}(x){f_{{Y_k}}}(y), \end{array} $

(10)

where $ {f_{X{Y_k}}} $ is the joint PDF of the kth distribution model, $ {c_{X{Y_k}}} $ is the copula density function for the kth component, and $ {f_{{X_k}}} $ and $ {f_{{Y_k}}} $ are the PDFs of X and Y for the kth component, respectively. The probability distribution of WS-SC can then be expressed as

$ {f_X}(x) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{f_{{X_k}}}(x)}, $

(11)

$ {F_X}(x) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{F_{{X_k}}}(x)}, $

(12)

and

$ {f_Y}(y) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{f_{{Y_k}}}(y)}, $

(13)

$ {F_Y}(y) = \sum\limits_{k = 1}^{{N_c}} {{w_k}{F_{{Y_k}}}(y)}, $

(14)

respectively.

To determine mixture bivariate statistical model parameters, the widely used EM algorithm proposed by Gupta and Chen (2011) was adopted. To derive the EM iteration function, a hidden random variable Z = {z_i}_{i = 1, 2, ···, N} (z_i ∈ {1, 2, ···, N_c}) is introduced, which indicates the particular bivariate distribution model from which the observations are sampled. Specifically, (x_i, y_i) is produced by $ {F_{X{Y_k}}} $ if z_i = k. Therefore, the joint PDF of (X, Y) and Z is

$ f\left[ {({x_i}, {y_i}), {z_i} = k|{\psi _k}} \right] = {w_k}{f_{X{Y_k}}}({x_i}, {y_i}|{\xi _k}) . $

(15)

For the kth component, the conditional PDF of Z can be expressed as

$\begin{align} f\left[ {{z_i} = k|({x_i}, {y_i}), {\psi _k}} \right] = &\frac{{f\left[ {({x_i}, {y_i}), {z_i} = k|{\psi _k}} \right]}}{{f\left[ {({x_i}, {y_i})|\psi } \right]}}\\ & = \frac{{{w_k}{f_{X{Y_k}}}({x_i}, {y_i}|{\xi _k})}}{{\sum\limits_{j = 1}^{{N_c}} {{w_j}{f_{X{Y_j}}}({x_i}, {y_i}|{\xi _j})} }} . \end{align} $

(16)

The iterant expression of the EM algorithm for the rth iteration can be derived from Eq. (16):

$\begin{array}{l} Q\left[ {\psi |{\psi ^{(r)}}} \right] = \sum\limits_{i = 1}^N {{E_{{z_i}|({x_i}, {y_i}), {\psi ^{(r)}}}}\left\{ {\ln f\left[ {({x_i}, {y_i}), {z_i} = k|{\psi _k}} \right]} \right\}}\\ = \sum\limits_{i = 1}^N {\sum\limits_{k = 1}^{{N_c}} {f\left[ {{z_i} = k\left| {({x_i}, {y_i})} \right., \psi _k^{(r)}} \right]\ln f\left[ {({x_i}, {y_i}), {z_i} = k\left| {{\psi _k}} \right.} \right]} }\\ = \sum\limits_{i = 1}^N {\sum\limits_{k = 1}^{{N_c}} {\frac{{w_k^{(r)}{f_{X{Y_k}}}\left({{x_i}, {y_i}\left| {\xi _k^{(r)}} \right.} \right)}}{{\sum\limits_{j = 1}^{{N_c}} {w_j^{(r)}{f_{X{Y_j}}}\left({{x_i}, {y_i}\left| {\xi _j^{(r)}} \right.} \right)} }}\ln \left[ {{w_k}{f_{X{Y_k}}}({x_i}, {y_i}\left| {{\xi _k}} \right.)} \right]} } . \end{array} $

(17)

The problem is then transformed into identification of ψ, which maximizes Eq. (17), to get the optimal estimation:

$\begin{array}{l} Q\left[ {\psi |{\psi ^{(r)}}} \right] = \sum\limits_{i = 1}^N {\sum\limits_{k = 1}^{{N_c}} {\frac{{w_k^{(r)}{f_{X{Y_k}}}\left({{x_i}, {y_i}\left| {\xi _k^{(r)}} \right.} \right)}}{{\sum\limits_{j = 1}^{{N_c}} {w_j^{(r)}{f_{X{Y_j}}}\left({{x_i}, {y_i}\left| {\xi _j^{(r)}} \right.} \right)} }}\ln {w_k}} } +\\ \;\;\;\;\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^{{N_c}} {\frac{{w_k^{(r)}{f_{X{Y_k}}}\left({{x_i}, {y_i}\left| {\xi _k^{(r)}} \right.} \right)}}{{\sum\limits_{j = 1}^{{N_c}} {w_j^{(r)}{f_{X{Y_j}}}\left({{x_i}, {y_i}\left| {\xi _j^{(r)}} \right.} \right)} }}\ln {f_{X{Y_k}}}\left({{x_i}, {y_i}\left| {{\xi _k}} \right.} \right)} }\\ \;\; = {Q_1}\left[ {w|{\psi ^{(r)}}} \right] + {Q_2}\left[ {\xi |{\psi ^{(r)}}} \right] .\end{array} $

(18)

The weight parameters w_k can be updated in agreement with Gibbs' inequality, which is given by

$ w_k^{(r + 1)} = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{w_k^{(r)}{f_{X{Y_k}}}\left({{x_i}, {y_i}\left| {\xi _k^{(r)}} \right.} \right)}}{{\sum\limits_{j = 1}^{{N_c}} {w_j^{(r)}{f_{X{Y_j}}}\left({{x_i}, {y_i}\left| {\xi _j^{(r)}} \right.} \right)} }}}, $

(19)

and the optimal value of ξ_k can be determined by setting the partial derivatives of $ {Q_2}[\xi |{\psi ^r}] $ to zero:

$ \frac{{\partial {Q_2}\left[ {\xi |{\psi ^{(r)}}} \right]}}{{\partial \xi }} = 0 . $

(20)

In practical applications, specifying suitable initial values for ψ is critical for the convergence of the EM algorithm. Random initial values are commonly used in such cases. The procedure involves randomly partitioning the data into N_c groups, representing the N_c components of the mixture model, followed by estimation of the parameters of the kth component using the subsample from the kth group. In this way, the initial values w_k can be expressed as

$ w_k^{(0)} = \frac{{{N_k}}}{N}, $

(21)

where N_k is the size of the kth subsample. Upon determining the initial values, the iterative EM algorithm can be performed by solving Eqs. (19) and (20) repeatedly. The iterative process is continued until the difference between two consecutive log-likelihood functions L(ψ) is less than a specified tolerance ε, demonstrating that the algorithm has converged and should be stopped:

$ \left| {L\left[ {{\psi ^{(r + 1)}}} \right] - L\left[ {{\psi ^{(r)}}} \right]} \right| < \varepsilon {\text{ (}}\varepsilon > 0{\text{)}}, $

(22)

where

$\begin{align} L(\psi) =& \sum\limits_{i = 1}^N {\log {f_{XY}}({x_i}, {y_i}|\psi)}\\ =& \sum\limits_{i = 1}^N {\log \left[ {\sum\limits_{k = 1}^{{N_c}} {{w_k}{f_{X{Y_k}}}({x_i}, {y_i}|{\xi _k})} } \right]} . \end{align}$

(23)

It should be worth noting that determining the magnitude of N_c is crucial for the application of the mixture bivariate statistical model. To determine the optimal N_c value, the BIC (Dziak et al., 2018) is used in this study

$ BIC = - 2L(\psi) + {N_\psi }\log (N), $

(24)

where N_ψ means the size of the model parameters. The optimal value of N_c can be determined by minimizing the BIC.

2.1.4 Error metrics

The efficacy of mixture distribution models was evaluated by applying the root-mean-square error (RMSE). Thus, the divergence between the joint probability densities of the original samples and simulated data from established models is measured by RMSE. To produce synthetic data, an algorithm (Liao et al., 2022; Yang et al., 2023) was utilized to generate bivariate random samples of a size greater than 10 times the number of observations. The bivariate space was then partitioned into M × L cells, where M and L refer to the number of bins in the WS-SC directions. Once the cells have been established, the number of bivariate observations (ED_ml) and simulated samples (SD_ml) that satisfy the constraints

$ \left\{ {\begin{array}{*{20}{l}} {\begin{array}{*{20}{c}} {{x_m} < x \leqslant {x_{m + 1}}, }&{m = 1, 2, \cdots, M} \end{array}} \\ {\begin{array}{*{20}{c}} {{y_l} < y \leqslant {y_{l + 1}}, }&{l = 1, 2, \cdots, L} \end{array}} \end{array}} \right., $

(25)

where (x_m, x_{m + 1}) and (y_l, y_{l + 1}) denote the boundaries of the cell with index (m, l), were computed for each cell. The RMSE for the bivariate distribution model can then be calculated using

$ RMSE = \sqrt {\frac{1}{{M \times L}}\sum\limits_{m = 1}^M {\sum\limits_{l = 1}^L {{{\left({\frac{{E{D_{ml}}}}{{{N_E}}} - \frac{{S{D_{ml}}}}{{{N_S}}}} \right)}^2}} } }, $

(26)

where N_E and N_S are the total number of measurements and simulated data, respectively. Generally, a smaller RMSE value indicates a superior fitting performance of the distribution model. It is worth noting, however, that the size of the divided cells may affect the RMSE values. In this work, the cell size was specified to be 0.5 m s⁻¹ × 0.005.

The fitting accuracy of the marginal distributions of WS-SC is proven by three metrics, namely, the squared Euclidean distance (D²), Kolmogorov-Smirnov distance (KS), and Pearson's correlation coefficient (R), which can be expressed as follows:

$ {D^2} = \frac{1}{N}\sum\limits_{i = 1}^N {{{(F_E^i - F_T^i)}^2}}, $

(27)

$ KS = \max \left| {{F_E} - {F_T}} \right|, $

(28)

$ R = \frac{{\sum\limits_{i = 1}^N {(F_E^i - \overline {{F_E}})(F_T^i - \overline {{F_T}})} }}{{\sqrt {\sum\limits_{i = 1}^N {{{(F_E^i - \overline {{F_E}})}^2}} } \sqrt {\sum\limits_{i = 1}^N {{{(F_T^i - \overline {{F_T}})}^2}} } }}, $

(29)

where N is the sample size, $ F_E^i $ is the empirical cumulative probability, $ F_T^i $ is the theoretical cumulative probability, and $ \overline {{F_E}} $ and $ \overline {{F_T}} $ are the means of $ F_E^i $ (i = 1, 2, ···, N) and $ F_T^i $ (i = 1, 2, ···, N), respectively. In general, smaller D² and KS values and higher R values mean better fitting performance.

2.2 Height-Dependent Wind Energy Assessment

Wind resource characteristic analysis was performed using the wind turbine model Aerodyne SCD 8.0/168 (ASCD-8). This model was selected to analyze the height-dependent wind resource characteristics. The basic characteristics of ASCD-8 are summarized in Table 1, and Fig.1 plots the power curve. Note that the hub height of this model is site-specific. Considering the prospective trend of increasing the wind turbine size in the future, hub heights of 60, 80, 100, 120, 140, 160, 180, and 200 m were specified for ASCD-8.

The WPD at height H_h is considered the main factor for wind resource characteristic analysis and is given by

$ WP{D_{{H_h}}} = \frac{1}{2}{\rho _{{H_h}}}V_{{H_h}}^3, $

(30)

where $ {\rho _{{H_h}}} $ is the air density at H_h and $ {V_{{H_h}}} $ is the wind speed at heights H_h. Considering the significant variation of air density with geographical location, the mean of air density from Yang et al. (2023) was used at different positions rather than the standard air density. To focus on the effect of wind shear on height-dependent wind resource characteristic analysis, air density is assumed to be constant with the change in hub height. Noteworthily, wind turbines are developed to produce electricity when wind speeds at hub height H_h are within V_cut-in and V_cut-out. The power law profile method was employed to determine the corresponding lower and upper limits of wind speed at H₀, denoted as V₁ and V₂, respectively:

$ {V_1} = \frac{{{V_{{\text{cut - in}}}}}}{{{{\left({{H_h}/H} \right)}^E}}}, $

(31)

$ {V_2} = \frac{{{V_{{\text{cut - out}}}}}}{{{{\left({{H_h}/H} \right)}^E}}} . $

(32)

Previous works (Liang et al., 2020; Wen et al., 2021; Yang et al., 2023) have modeled wind speeds at various hub heights individually to determine the height-dependent WPD. However, in this study, the improved WPD at any arbitrary H_h can be acquired using the established mixture bivariate statistical models:

$\begin{array}{l} WP{D_{{H_h}}} = \\ \;\;\frac{1}{2}\int_0^\infty {\int_{{V_1}}^{{V_2}} {{\rho _{{H_h}}}{{\left[ {{V_{{H_0}}}{{\left({\frac{{{H_h}}}{{{H_0}}}} \right)}^E}} \right]}^3}{f_{XY}}({V_{{H_0}}}, E)\text{d}{V_{{H_0}}}} } \text{d}E . \end{array} $

(33)

EP is another important factor in wind resource characteristic analysis (Lin and Dong, 2016). A higher magnitude of EP means that more wind electricity can be produced. In terms of offshore wind power projects, a precise estimation of EP is needed to ensure economic feasibility:

$ EP = T \times {P_w}({V_{{H_h}}}), $

(34)

where T is the number of hours over a specified period and P_w(·) is the power curve shown in Fig.1. Using the developed mixture bivariate distribution model, Eq. (34) can be rewritten as

$ EP = T\int_0^\infty {\int_{{V_1}}^{{V_2}} {{P_w}({V_{{H_h}}}){f_{XY}}({V_{{H_0}}}, E)} } {\text{ d}}{V_{{H_0}}}{\text{d}}E . $

(35)

The power curves are applicable under standard air density conditions. However, in practical applications, air density affects the wind resource evaluation and should be considered (Jung and Schindler, 2019a; Liang et al., 2020; Huang et al., 2023; Yang et al., 2023). Thus, to estimate wind speed under standard air density (V_st), it is imperative to transform wind speed under real air density $ ({V_{{H_h}}}) $ using the Svenningsen approach (Svenningsen, 2010)

$ {V_{st}} = {V_{{H_h}}}{\left({\frac{{{\rho _{{H_h}}}}}{{{\rho _{st}}}}} \right)^{z({V_{{H_h}}})}}, $

(36)

where ρ_st = 1.225 kg m⁻³ is the standard air density and $ z({V_{{H_h}}}) $ is a correction variable:

$ z({V_{{H_h}}}) = \left\{ {\begin{array}{*{20}{c}} {\frac{1}{3}}&, &{{V_{{H_h}}} \leqslant 8} \\ {\frac{1}{3} + \frac{1}{3}\frac{{{V_{{H_h}}} - 8}}{5}}&, &{8 < {V_{{H_h}}} < 13} \\ {\frac{2}{3}}&, &{{V_{{H_h}}} \geqslant 13} \end{array}} \right. . $

(37)

Then, V_st can be computed by combining Eq. (36) and the power law profile method equation (Eq. (39)) (Gualtieri, 2019):

$ {V_{st}} = {V_{{H_0}}}{\left({\frac{{{H_h}}}{{{H_0}}}} \right)^E}{\left({\frac{{{\rho _{{H_h}}}}}{{{\rho _{st}}}}} \right)^{z({V_{{H_h}}})}}, $

(38)

$ {V_{{H_h}}} = {V_{{H_0}}}{\left({\frac{{{H_h}}}{{{H_0}}}} \right)^E}, $

(39)

where $ {V_{{H_h}}} $ and $ {V_{{H_0}}} $ are the wind speeds at heights H_h and H₀, respectively, and E is the wind shear coefficient. Then, EP can be computed using the following expression:

$ EP = T\int_{{E_{\min }}}^{{E_{\max }}} {\int_{{V_1}}^{{V_2}} {{P_w}\left[ {{V_{{H_0}}}{{\left({\frac{{{H_h}}}{{{H_0}}}} \right)}^E}{{\left({\frac{{{\rho _{{H_h}}}}}{{{\rho _{st}}}}} \right)}^{z({V_{{H_h}}})}}} \right]{f_{XY}}({V_{{H_0}}}, E)} } {\text{ d}}{V_{{H_0}}}{\text{d}}E . $

(40)

To calculate the annual energy production (AEP), T is defined as 365 × 24 = 8760. Meanwhile, to compute the seasonal energy production (SEP), T is defined as 8760/4 = 2190.

Moreover, in practical wind energy applications, CF is a crucial technical measure of the capacity of wind turbines to convert wind resources into electricity and is defined as the ratio of the AEP of a wind turbine to the rated power output over the span of one year (Quan and Leephakpreeda, 2015). Mathematically,

$ CF = \frac{{AEP}}{{T \times {P_R}}}, $

(41)

where T = 8760.

2.3 Estimation of the Influence of Wind Shear Variation

To underscore the influence of wind shear coefficient variability on height-dependent wind resource characteristics, the percentage change in wind resource assessment factors (PC) was introduced and computed using

$ PC = \frac{{{\eta _c} - \eta }}{\eta } \times 100\%, $

(42)

where η and η_c are the wind resource assessment factors for real wind shear coefficient conditions and the assumed constant wind shear coefficient condition (E = 0.11), respectively. Using the resulting marginal distribution of wind speed, the WPD at height H_h for constant wind shear conditions was calculated using

$ WP{D_c} = \frac{1}{2}\int_{{V_1}}^{{V_2}} {{\rho _{{H_h}}}{{\left[ {{V_{{H_0}}}{{\left({\frac{{{H_h}}}{{{H_0}}}} \right)}^E}} \right]}^3}{f_X}({V_{{H_0}}})} {\text{ }}\text{d}{V_{{H_0}}} . $

(43)

Similarly, EP and CF for a constant wind shear condition can be derived from

$ E{P_c} = T\int_{{V_1}}^{{V_2}} {{P_w}\left[ {{V_{{H_0}}}{{\left({\frac{{{H_h}}}{{{H_0}}}} \right)}^E}{{\left({\frac{{{\rho _{{H_h}}}}}{{{\rho _{st}}}}} \right)}^{z({V_{{H_h}}})}}} \right]{f_X}({V_{{H_0}}}){\text{ }}} {\text{d}}{V_{{H_0}}}, $

(44)

$ C{F_c} = \frac{{AE{P_c}}}{{T \times {P_R}}} . $

(45)

Let PC_WPD, PC_EP, and PC_CF be the percentage changes in WPD, EP, and CF, respectively. A positive PC value refers to the overestimation of these factors due to the assumption of a constant wind shear coefficient, and a negative value means underestimation.

3 Data Description and Analysis

To examine the adaptability of the designed bivariate statistical model to represent diverse wind regimes and to study the spatial characteristics of offshore wind resources, research was conducted using samples taken from six offshore wind power projects in China to represent different wind regimes in different sea areas; detailed information is presented by Yang et al. (2023). Fig.2 displays a map of the geographical distribution of the selected sites.

The realization of a reliable and accurate wind energy evaluation can be improved by using a longer time series with high-quality data. In our work, a continuous 30-year hourly reanalysis dataset during 1993 – 2022, with the horizontal 10 m wind speed and wind shear coefficient at six positions, was taken from ERA5, which is the most contemporary reanalysis produced by the ECMWF that has undergone validation and finds widespread application in wind energy research (Soukissian et al., 2021; Patel et al., 2022; Wu et al., 2022b). The spatiotemporal resolutions for the wind parameters are 0.25˚ × 0.25˚ and 1 h (Rivas and Stoffelen, 2019; Hersbach et al., 2020).

The horizontal 10 m wind speed (V₁₀) was calculated using

$ {V_{10}} = \sqrt {u_{10}^2 + v_{10}^2}, $

(46)

where u₁₀ and v₁₀ are the 10 m u-component of wind and the 10 m v-component of wind, respectively, obtained from the single level product of ERA5.

The wind shear coefficient cannot be directly extracted from ERA5. The literatures (Jung and Schindler, 2017; Han et al., 2022; Debnath et al., 2023) suggested that the logarithm of both sides of the power law profile method equation (Eq. (39)) can be taken to compute the wind shear coefficient

$ E = \frac{{\log ({V_{{H_h}}}) - \log ({V_{{H_0}}})}}{{\log ({H_h}) - \log ({H_0})}} . $

(47)

In this method, the wind shear coefficient can be estimated on the basis of wind speed data at two different heights and Eq. (47).

He et al. (2022) demonstrated that the power law profile method can be employed at heights less than 400 m in the atmospheric boundary layer (ABL). To acquire a more robust estimation of the wind shear coefficient, all available wind speed data with heights not greater than 400 m in the ABL were adopted in our work. The detailed procedure is presented below.

Step 1. Calculate 100 m wind speed (V₁₀₀) based on the 100 m u-component and v-component of wind taken from ERA5.

Step 2. Calculate wind speeds at different pressure levels (p_k (k = 1, 2, 3, 4) = 925, 950, 975, and 1000 hPa) denoted as V_k (k = 1, 2, 3, 4) using the corresponding u-component and v-component of wind taken from the pressure level products of ERA5.

Step 3. Determine the heights H_k (k = 1, 2, 3, 4) of V_k (k = 1, 2, 3, 4) using Eq. (48), which is derived from the barometric law given in Eq. (49) (Jung and Schindler, 2019a).

$ {H_k} = \left[ {1 - {{\left({\frac{{{p_k}}}{{{p_{{H_0}}}}}} \right)}^{\frac{1}{{5.2561}}}}} \right] \times \frac{{{T_{{H_0}}}}}{{0.0065}}, {\text{ }}k = 1, 2, 3, 4, $

(48)

$ {p_{{H_h}}} = {p_{{H_0}}}{\left({1 - \frac{{0.0065 \cdot \Delta h}}{{{T_{{H_0}}}}}} \right)^{5.2561}}, $

(49)

where $ {p_{{H_0}}} $ and $ {T_{{H_0}}} $ are the sea surface pressure and temperature, respectively.

Step 4. At each hour of the investigated period, select hourly wind speed data, the heights of which are less than the ABL height (obtained from ERA5) and 400 m, from V₁₀, V₁₀₀, and V_k (k = 1, 2, 3, 4). Let N_w be the number of selected wind speed samples. If N_w is greater than 2, the least squares fitting method is employed using Eq. (47) for wind shear coefficient estimation. Otherwise, calculate this coefficient using two wind speeds with lower heights.

Boxplots of the total WS-SC samples at the six positions are depicted in Fig.3. Both WS-SCs show significant differences in their probabilistic distribution characteristics across various geographical positions. These variations imply that the wind resources at the six positions may differ. Specifically, position L4 records the highest wind speed median, whereas position L2 has the lowest. Position L1 has the highest wind shear coefficient median, whereas position L5 has the lowest. Note that the majority of wind shear coefficients at all chosen positions are lower than the assumed constant value (the black line in Fig.3) for offshore conditions due to different complex mechanisms (e.g., wind-wave interaction). Thus, statistical modeling of the wind shear coefficient is necessary for accurate height-dependent wind resource characteristic analysis.

Previous research revealed that wind energy exhibits significant seasonal variation (Shu et al., 2015; Yang et al., 2023), and the wind shear coefficient also demonstrates strong seasonality (Yan et al., 2022). Therefore, it is imperative to perform a seasonal investigation of height-dependent wind resource characteristics. The boxplots of the seasonal WS-SC samples are illustrated in Fig.4, where both variables show strong seasonal dependence. The highest wind speed median is observed in winter for all locations, whereas the lowest value is observed in spring or summer. The seasonal pattern of the wind shear coefficient is different, with the peak value in spring or summer and the minimum in autumn or winter. Season-dependent distribution characteristics emphasize the importance of seasonal analysis of wind resources.

4 Results and Discussion 4.1 Mixture Bivariate Statistical Model of Wind Speed and Wind Shear

A mixture bivariate statistical model was developed in this work to capture the multimodal structures exhibited in the joint probability distribution of WS-SC. The number of components of the models for different samples was identified based on BIC. The model parameters were computed using the EM algorithm. To assess the performance of the proposed model, bivariate samples were generated. The RMSE was computed to measure the differences between the empirical data and the simulated bivariate samples, and the results are summarized in Table 2. The joint probability density plots of the WS-SC for the total samples and simulated datasets are presented in Figs.5 and 6, respectively. Similarly, Figs.7 and 8 show the joint probability density plots for seasonal data.

The marginal distributions of WS-SC play a critical role in practical offshore wind energy industry applications. Upon establishing the mixture bivariate statistical models, the corresponding marginal distributions for WS-SC can be derived from Eqs. (11) and (13). The marginal distributions of WS-SC for the total and seasonal samples are shown in Figs.9 − 12. To compare with the developed model, histograms of the empirical probability density distributions are also included.

The lower RMSE values exhibit the remarkable fitting performance of the mixture bivariate model for the joint probability density distribution of WS-SC in both the total and seasonal samples. As shown in Figs.5 and 6, the peaks of the modes in the wind shear coefficient direction are widely separated, demonstrating distinct multimodal characteristics in the marginal distribution of the wind shear coefficient, as displayed in Fig.10. In contrast, there is only one hump in the probability density distributions of wind speed, as shown in Fig.9. This phenomenon occurs because the components in this direction are close to each other. The bivariate joint probability density plots of the empirical and simulated samples disclose that the developed mixture distribution model accurately captures the multimodal nature of the joint probability distribution of WS-SC.

In agreement with the distribution characteristics of WS-SC displayed in Fig.4, the probabilistic structures of these two variables show conspicuous seasonality. This phenomenon is well displayed in Figs.7 and 8. The modes of the joint distributions in the wind shear coefficient direction during autumn and winter are sufficiently near each other, resulting in unimodal marginal distributions. In contrast, the constituent components of the wind shear coefficient are far from each other, leading to the existence of two humps in the marginal distributions of this variable. The developed mixture models for seasonal data can offer adequate fitting performance for the joint distributions, further ensuring the accuracy of the seasonal wind energy characteristic analysis.

The suitability of the mixture Weibull and mixture normal distributions for WS-SC is shown in Figs.8 – 11. The marginal distributions of the wind shear coefficient for both the total and seasonal samples show multimodal characteristics that cannot be described by a single parametric model. Therefore, a linear combination of conventional distributions is needed (Huang and Dong, 2019), and the mixture normal distribution proves to be an effective fit. Note that although the probability density distributions of wind speed typically display unimodal shapes, their primary characteristics remain satisfactorily described by the mixture Weibull distribution. Thus, irrespective of whether the empirical probability density distributions are unimodal or multimodal, the mixture models can offer a suitable fit. The test statistics for the total WS-SC samples are summarized in Table 3. It can be easily observed that the established mixture univariate distributions are capable of providing outstanding fitting performances. The mixture univariate distributions for the seasonal data also show satisfactory performance, based on the results in Table 4.

4.2 Height-Dependent Wind Resource Characteristic Analysis

Figs.13 and 14 illustrate the height-dependent annual average WPD and seasonal average WPD, respectively. It is observed that WPD increases with increasing hub height at the investigated positions. Specifically, at L4, the annual WPD variation is greater than 200 W m⁻² with changing hub height from 60 to 200 m, whereas at L6, this variation is approximately 100 W m⁻². These different effects are associated with the joint distribution characteristics of WS-SC. Thus, an accurate bivariate statistical model is imperative for an accurate assessment of height-dependent WPD. The seasonal variations in WPD resulting from the changes in hub height are also significant. At L1, the seasonal WPD variation for spring can be up to 400 W m⁻², whereas that for autumn and winter is only about 100 W m⁻² at L1, L2, L3, and L6. These differences strongly underscore that performing a seasonal investigation of wind resource characteristics is crucial. In addition, influenced by both WS-SC, the WPD at a specific height varies significantly with space. For example, L4 and L5 have richer wind resources with annually average WPD exceeding 600 at the height of 60 m. In contrast, L1, L2, L3, and L6 have a poorer potential for offshore wind resource extraction, with an annual average WPD of less than 600 W m⁻² even at the height of 200 m. As for seasonal average WPD, L1 and L2 have greater potential for wind energy in spring, whereas winter is the period with particularly abundant wind resources at L3, L4, L5, and L6.

Figs.15 and 16 illustrate the height-dependent AEP and SEP of the selected wind turbine model. The results disclose a positive correlation between hub height and wind electricity generation, which is consistent with the WPD analysis. For instance, at L2, L3, and L4, increasing the hub height from 60 to 200 m can add approximately 5 GW h to the AEP. Similarly, peak values of variation in SEP are found in spring at L1 and L2, amounting to over 2 GW h. Notably, the rate of change in EP decreases with increasing hub height, gaining critical insights into whether the additional EP produced by a higher hub height can offset the incremental cost related to wind turbine design, manufacture, installation, and other factors, ultimately allowing cost-efficient decision-making (Jung and Schindler, 2017). The significance of this analysis is thus emphasized. Moreover, the results reveal distinct spatial and seasonal variations in the height-dependent EP. Specifically, L4 records the highest AEP, whereas L2 records the lowest. The wind turbine generates the most energy in winter at L3, L4, L5, and L6, whereas the highest SEP values are found in spring at L1 and L2. Meanwhile, the EP during summer is generally lower along the eastern coast of China.

From Eq. (41), the vertical variation and spatial distribution of the CF would coincide with those of the AEP, as presented in Table 5. Particularly, the CF increases in proportion to the hub height, with more fluctuations taking place at L2, L3, and L4, where the variation in CF results in approximately 0.08 with increasing hub height from 60 to 200 m. It is worth noting that the CF for a hub height of 200 m at L1, L2, L3, and L6 is still below 0.45, showing that the wind turbine is particularly well suited for electricity generation in these areas (de Sousa Gomes et al., 2019). Also, coinciding with the spatial distribution of the height-dependent AEP, the peak and valley values of the CF for the wind turbine are determined at L4 and L2, respectively.

4.3 Influence of Wind Shear on the Wind Resource Assessment

Table 6 lists the magnitudes of PC_WPD for the annual average WPD, and Fig.17 presents those for the seasonal average WPD. It is shown that PC_WPD generally increases with the increase in hub height for both the annual and seasonal average WPD. The median of the wind shear coefficient in Fig.3 is less than the assumed constant value, leading to negative PC_WPD values, depicting an overestimation of the annual averaged WPD at all six positions. The maximum PC_WPD locates at L5 with a magnitude of 52.93% at a height of 200 m due to the largest deviation between the real and assumed wind shear coefficients. L1 shows the least PC_WPD, with values less than 7.5% at all studied hub heights. In addition, impacted by the seasonal variation of the wind shear coefficient, the maximum PC_WPD for the seasonal average WPD is concentrated on autumn and winter. Except for the WPD of spring and summer at L1 and L2, the seasonal WPD is usually overestimated at the considered locations.

Table 7 and Fig.18 exhibit the results of PC_EP for AEP and SEP, respectively. It is clear that the wind shear coefficient plays a significant role in the calculation of AEP and SEP, in particular with increasing hub heights in the future. For a particular hub height, the maximum PC_EP values for the total data are observed at L6, followed by L3 or L5, whereas the minimum value is observed at L1 or L4. The PC_EP values for the seasonal data are generally higher in autumn and lower in summer. Notably, a significant PC_EP can seriously influence the economic feasibility assessment of offshore wind power projects in practical applications. Moreover, similar to the seasonal data results of PC_WPD, the variation of the wind shear coefficient also influences SEP differently in the northern and southern regions. For example, at L1, the SEP values of spring and summer are underestimated when assuming a constant wind shear coefficient, whereas at L3, L4, L5, and L6, the SEP values of all seasons are overestimated.

In addition, based on Eqs. (44) and (47), it can be inferred that the PC_CF results are identical to those of PC_EP for the total data. Hence, the characteristics of the wind shear coefficient at the specified sites must be thoroughly investigated for accurate evaluation of the efficiency of the wind turbine and to determine if the selected wind turbine is a good fit for local wind conditions (Koholé et al., 2023).

5 Conclusions

In this work, we investigated the height-dependent wind resource characteristics at six positions in the coastal sea area of China, where offshore wind farms are being planned for development. A mixture bivariate statistical model was built to model the joint distribution of WS-SC under mixed atmospheric conditions. For each component, the marginal distributions of the WS-SC were fitted using a Weibull distribution and a normal distribution, respectively, while the dependency structure was described by a Gaussian copula. The optimal number of mixture model components was determined using the BIC, and the model parameters were computed on the basis of the EM algorithm. The developed model was corroborated using both total and seasonal reanalysis data at the considered positions.

The joint probability density distributions of the WS-SC for both the total and seasonal samples at the investigated coastal positions present distinct multimodal characteristics. The developed joint statistical model can adequately capture these complex statistical patterns. Furthermore, the marginal distributions of WS-SC, whether unimodal or multimodal, can be fitted well by the proposed mixture of Weibull and normal distributions. Thus, it is suggested that mixture bivariate and univariate statistical models be utilized to model wind parameters and comprehensively evaluate wind resource potential. Also, the proposed model is an effective approach to analyzing wind resource characteristics at hub heights of 10 – 200 m. It enables the evaluation of wind resource characteristics at any arbitrary hub height without requiring measured data at a specific hub height.

The proposed model was used to assess the height-dependent wind resource characteristics at investigated positions. WPD, EP, and CF were evaluated at different hub heights, and the results demonstrate that these wind energy assessment factors increase with increasing hub height across all positions. Notably, spatial and temporal variability is found in these factors, showcasing the significance of a comprehensive spatial and seasonal investigation of height-dependent wind resources.

Moreover, the influence of wind shear coefficient variation on wind resource characteristics was quantitatively analyzed. The results demonstrated that the errors resulting from assuming a constant wind shear coefficient when calculating WPD, EP, and CF cannot be ignored, especially with the increase in hub height in the future. It was observed that these factors were overestimated at the studied positions, and the percentage change of these factors varied between the northern and southern regions. Also, the percentage change was generally higher in autumn and lower in summer for seasonal factors. Provided the potential consequences of misestimating wind resource characteristics on the economic feasibility of offshore wind power projects, proper consideration of wind shear coefficient variation and an accurate statistical model are imperative in wind energy research.

However, this work has some limitations. First, because of the lack of adequate measurements of wind speed at different heights, the wind shear coefficient was estimated using reanalysis wind speeds according to the developed procedure. In addition, only reanalysis wind speed and wind shear coefficient data at six locations along the Chinese coastline were adopted to validate the developed mixture bivariate distribution model without considering its applicability to distinct sea areas of the world. In the future, we will collect long-term wind speed observations at different heights to compute the wind shear coefficient and correct the values estimated by the reanalysis data. In addition, the general applicability of the proposed model will be shown well using more data distributed worldwide, and the influence resulting from the spatial variation of the wind shear coefficient on wind energy assessment will be further studied.

Acknowledgements

This study was supported by the Key R & D Program of Shandong Province, China (No. 2021ZLGX04), and the National Natural Science Foundation of China (No. 52171284). The authors would like to acknowledge the ECMWF for providing the ERA5 dataset, which is available at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5. The authors would like to acknowledge the website http://en.wind-turbine-models.com/ for providing the adopted wind turbine model.

Nomenclature

AEP, annual energy production

BIC, Bayesian information criterion

C [∙], copular function

c [∙], copula density function

CF, capacity factor

D², squared Euclidean distance

E, wind shear coefficient

EP, energy production

F_XY(·), bivariate CDF of X and Y

f_XY(·), bivariate PDF of X and Y

F_X(·), CDF of X

f_X(·), PDF of X

F_Y(·), CDF of Y

f_Y(·), PDF of Y

H_h, hub height

KS, Kolmogorov-Smirnov distance

L(·), log-likelihood function

N, sample size

N_c, number of components in the mixture model

N_ψ, size of the parameters

PC, percentage change

PC_CF, percentage change in CF

PC_EP, percentage change in EP

PC_WPD, percentage change in WPD

$ {p_{{H_h}}} $, air pressure at H_h

P_R, rated power of the wind turbine

P_w(·), power output

R, Pearson's correlation coefficient

RMSE, root-mean-square error

SEP, seasonal energy production

T, number of hours over a specified period

$ {T_{{H_h}}} $, temperature at H_h

$ {T_{{H_0}}} $, temperature at height H₀

V, wind speed

V_cut-in, cut-in wind speed

V_cut-out, cut-out wind speed

$ {V_{{H_h}}} $, wind speed at hub height H_h

V_st, wind speed for the standard air density condition

w, weight parameter

WPD, wind power density

x, wind speed

y, wind shear coefficient

z(·), correction variable

α, shape parameter

β, scale parameter

η, mean of y

σ, standard deviation of y

$ {\rho _{{H_h}}} $, air density at H_h

ρ_st, standard air density (1.225 kg m⁻³)

Φ(·), CDF of the standard normal distribution

γ, parameter of the Gaussian copula function

Author Contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Zihao Yang. The draft of the manuscript was written by Zihao Yang, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Data Availability

All data generated and analyzed during this study are included in this published article and its additional files.

Declarations

Ethics Approval and Consent to Participate

This article does not contain any studies with human participants or animals performed by any of the authors.

Consent for Publication

Informed consent for publication was obtained from all participants.

Conflict of Interests

The authors declare that they have no conflict of interests.