Joint Modelling of Wave Energy Flux and Wave Direction

: In the context of wave resource assessment, the description of wave climate is usually conﬁned to signiﬁcant wave height and energy period. However, the accurate joint description of both linear and directional wave energy characteristics is essential for the proper and detailed optimization of wave energy converters. In this work, the joint probabilistic description of wave energy ﬂux and wave direction is performed and evaluated. Parametric univariate models are implemented for the description of wave energy ﬂux and wave direction. For wave energy ﬂux, conventional, and mixture distributions are examined while for wave direction proven and efﬁcient ﬁnite mixtures of von Mises distributions are used. The bivariate modelling is based on the implementation of the Johnson–Wehrly model. The examined models are applied on long-term measured wave data at three offshore locations in Greece and hindcast numerical wave model data at three locations in the western Mediterranean, the North Sea, and the North Atlantic Ocean. A global criterion that combines ﬁve individual goodness-of-ﬁt criteria into a single expression is used to evaluate the performance of bivariate models. From the optimum bivariate model, the expected wave energy ﬂux as function of wave direction and the distribution of wave energy ﬂux for the mean and most probable wave directions are also obtained.


Introduction
Marine renewable energy applications (especially, offshore wind, wave, and current/tide) are gaining ground as means for deriving electrical power due to the commitment of the EU member states to reduce pollutant emissions by substituting fossil fuels with clean energy sources [1,2]. To date, EU funds have supported particularly the wave energy sector with respect to technology R&D. Although the technological progress is remarkable with numerous design concepts and prototypes to harness energy from sea waves, the corresponding commercial stage follows a slower pace.
The particular advantages of wave energy (e.g., high energy density, small energy loss, predictability) render its exploitation promising not only in areas with high wave energy flux like the Atlantic Ocean but also in areas with relatively calm wave climate and thus with intermediate levels of power availability, like the Mediterranean Sea; see, for example, [1,3,4]. Usually, the later areas are characterized by low variability achieving better survivability, without increase of capital expenditure [5]. In order to make wave energy production feasible and economically viable, it is necessary, among others, to deploy wave energy converters (WECs) based on the accurate knowledge of the wave climate and its temporal variability in the area of interest [6]. The wave climate variability is also responsible to a great extent for the uncertainty in electricity production of a wave energy project, an aspect that has been studied by [7,8].
Such studies are usually confined in the statistical description of particular wave spectral parameters, such as significant wave height H S defined as H S ∼ = 4 √ m 0 ; energy period T e , T e = m −1 /m 0 ; and mean wave direction Θ W = tan −1 S ηη ( f ,θ) sin θdθd f S ηη ( f ,θ) cos θdθd f , where m n = ∞ 0 f n S ηη ( f )d f , n = · · · , −2, −1, 0, 1, · · · , is the n-order spectral moment, S ηη ( f ) is the spectral density function, and f is the spectral frequency; see, e.g., [15]. However, there are recent studies on wave energy resource characterization [16] that are based on the six IEC parameters, which include at least two directional parameters: direction of maximum energy and directionality coefficient.
Wave energy flux per unit length of wave front (W/m), which is also known as wave power density or wave energy potential, can be evaluated from the following relation: where ρ is the density of seawater, g is the acceleration of gravity, c g ( f , d) is the group velocity, d is the water depth, and θ is the direction of wave propagation; see also [17,18]. For deep water, the above relation for the total wave energy flux is simplified as follows: which is usually expressed through H S and T e via the following relation: where ρ = 1025 kg/m 3 and g = 9.8066 m/s 2 .
Regarding wave climate analysis per se, a lot of work has been done for the modelling of H S and T e (or spectral peak period T P and mean zero-up crossing period T Z ). Among the most common probability distributions for modelling H S (and/or T e /T P /T Z ) are the Beta, the Gamma, the Generalized Gamma, the Weibull, and the Lognormal distributions, with the latter one being the most widely used. An alternative modelling approach is non-parametric methods; in this context, kernel density estimation is frequently applied. Essentially, the kernel is used as a weighting function centred at the data points and its extension (around the data points) is defined by a smoothing parameter. Some standard reference books for non-parametric estimation are [25,26].
However, the univariate models are not sufficient in the context of a detailed longterm analysis, especially with respect to design purposes, load calculations, and reliability analyses; thus, the bivariate description of the above-mentioned wave parameters is necessary; [27][28][29]. Different approaches have been presented for the joint description of wave height and period. For instance, parametric bivariate distributions, conditional-distribution approach, kernel-based models, and copula-based models as well, provide alternative probabilistic description of these variables; see, e.g., [27,[30][31][32]. Furthermore, [33,34] have provided distributions for the joint description of wave power-wave height and wave power-wave period for individual waves, an approach that is complementary in long-term analysis of the wave climate at a particular site.
WECs can be classified regarding their direction of orientation (and size) into three categories: attenuators, point absorbers, and terminators; see, e.g., [35]. Unlike a point absorber that converts energy from waves coming from all directions, attenuators and terminators should be elongated parallel and perpendicular to the propagation of wave direction, respectively, to yield high-energy conversion efficiency. This sensitivity to directionality affects the standard deviation of the capture length (i.e., ratio of the WEC Processes 2021, 9, 460 3 of 24 power output to the sea state power) of WECs as has been demonstrated by [8]. Moreover, as is emphasized in [36], apart from H S and T e , the performance of a WEC should be based on additional attributes of the sea state including variability of mean wave direction, directional spread and spectral shape, and that the uncertainties of the performance results of the WECs can be reduced if these dependencies are explored and quantified; see also [37]. In the same context, [38] highlighted that developers of WECs should provide information regarding the effects of wave direction on the performance of WECs for an accurate assessment of wave energy since it affects wave propagation and in turn, the distribution of wave resource. A similar suggestion has been provided earlier by [36], advising to provide power matrix for different mean wave directions with the related uncertainty matrix. Reference [7] also has noted that vicinal sites characterized by similar resource levels may exhibit differences in acquirable wave energy with time due to the varying wave directional exposure.
The importance of analysing wave direction for performance optimization purposes of a WEC, and generally, in wave energy studies, has been also emphasized in [39,40]. In [41], variability of wave heights and wave power has been assessed with respect to wave direction for a more realistic wave resource analysis at different locations in the southwest UK. The directional dependence of wave height and, consequently, wave energy generation was more than evident. The authors concluded that site selection for commercial development of wave energy, reduction of grid integration hindrances, and increase of power smoothing are highly dependent on directional exposure of sites. Reference [42] has assessed wave power and its temporal variability for various water depths off the central west coast of India and examined the directional distribution of wave power at shallower water depths. Reference [43] has highlighted the wave (and wind) directionality effects on offshore wind turbines, supported by a jacket structure during hurricane events, while [44,45] have stressed the significance of wave direction, among others, that may lead to larger fatigue damage for certain platforms under specific conditions. Overall, the accurate description of wave directional characteristics and their variability is important when assessing wave energy potential at a site; see, e.g., [46], while wave direction affects significantly wave propagation, and in turn, wave resource distribution, in an area.
In this work, instead of modelling the involved individual random variables in wave energy assessment studies, i.e., significant wave height H S , energy period T e , and mean wave direction Θ W (Θ W is subsequently denoted as Θ), the joint assessment of "wave energy flux P" and "mean wave direction Θ" is performed. To the authors' knowledge, the joint description of P and Θ has not been presented yet, although it is of high importance for the emerging wave energy sector and relevant wave resource analysis applications. Moreover, some straightforward applications derived from the proposed methodology are also presented.
In this framework, numerous univariate conventional and mixture parametric models have been addressed with respect to the (linear) variable P. The directional variable Θ is adequately described through a finite mixture of the von Mises distribution that has already proven to be efficient in this case [40,47,48]. The construction of the bivariate distribution functions for (P, Θ) is accomplished through the Johnson-Wehrly (JW) parametric model. An important feature of the JW model is that its density functions rely on the corresponding univariate marginal distributions while an additional parameter quantifies the correlation/dependence of the linear and directional variable. A detailed evaluation of the resulting bivariate distributions is made by applying five bin-defined metrics, namely, adjusted coefficient of determination, index of agreement, mean absolute error, root mean square error, and χ 2 statistic, as well as a synthetic criterion that combines the above criteria into a single expression. See also [47,49,50].
The structure of the present work is the following: in Section 2, wave data obtained from three buoy locations in the Greek Seas along with data from numerical wave models in three locations (Gulf of Lion, Gulf of Biscay, and North Sea) are described, and a preliminary statistical analysis is performed for wave energy flux P and mean wave direction Θ (subsequently referred to as simply wave direction). The examined univariate probability models for P and Θ are presented in Section 3. Specifically, for the parametric modelling of wave energy flux, several conventional and mixture distributions are applied while wave direction is modelled by using a finite mixture of the von Mises distribution. The bivariate JW probability model is briefly presented in Section 4. Section 5 deals with the statistical criteria used to evaluate the performance of both univariate and bivariate models. Based on the goodness-of-fit evaluation criteria, the derived bivariate models are compared, and the corresponding results are discussed in Section 6. In Section 7, the best bivariate model for (P, Θ) for each examined area is implemented to provide analytic results that could be used in feasibility studies regarding wave farm development in the corresponding areas. The last section includes the concluding remarks of this analysis.

Description of Wave Data
Three different sources of wave data have been used in this work, corresponding in offshore areas with different wave characteristics and in turn, wave climate. The first source is in situ wave data from the Greek Seas (Aegean and Ionian seas), obtained from three oceanographic buoys of the POSEIDON marine monitoring network [51,52]. Wave measurements from buoys are confined in the Greek Seas as this data source was available to the authors. The wave parameters that were obtained were the significant wave height H S and the energy period T e in order to estimate wave energy flux P, and the wave direction Θ. The geographical coordinates, the measurement period, the sample size, and the water depth for each buoy location are listed in Table 1; see also Figure 1. The wave data were firstly filtered while only concurrent wave measurements of the three parameters (H S , T e , and Θ W ) were taken into consideration. Furthermore, wave data obtained from numerical wave models were also analysed. The first numerical wave model data set is produced by the Mediterranean Sea Waves forecasting system, i.e., a wave model that is based on WAM Cycle 4.5.4 and is provided by the EU Copernicus service [53] with spatial resolution 0.042 • × 0.042 • . The selected grid points correspond to locations in the gulfs of Lion (France) and Biscay (Spain). The second wave data set is obtained from the ERA20C reanalysis product of the European Centre for Medium-Range Weather Forecasts (ECMWF), available at 1.25 • horizontal spatial resolution. Ocean waves are described by spectra on 25 frequencies and 12 directions [54]. The selected grid point corresponds to a location offshore southern Norway, in the North Sea. All the spectral parameter time series (defined in Section 1) that are analysed in the present work have a 3-h step. See also Table 1 and Figure 1. More details and discussion as regards spatial distribution of wave spectral parameters, numerical wave models, and their validation through field data can be found for example in [55][56][57].  The following linear and directional parameters are necessary for the computations regarding the probability density functions presented in Sections 3 and 4. Specifically, the statistical parameters of wave energy flux , i.e., mean , median 0.5 and maximum values , standard deviation , and coefficient of variation = ⁄ , where is the mean value of , skewness , and kurtosis are presented in Table 2. For the measured data, the highest mean, median, maximum, standard deviation, and skewness values are recorded at Mykonos, and for the model wave data, the highest values for median, maximum, and standard deviation are encountered at Gulf of Biscay. Gulf of Lion exhibits by far the highest , skewness, and kurtosis values. Clearly, the most energetic wave climate from the examined locations is encountered in the Gulf of Biscay.  The following linear and directional parameters are necessary for the computations regarding the probability density functions presented in Sections 3 and 4. Specifically, the statistical parameters of wave energy flux P, i.e., mean µ P , median P 0.5 and maximum values max P , standard deviation σ P , and coefficient of variation CV P = σ P /µ P , where µ P is the mean value of P, skewness Sk, and kurtosis Ku are presented in Table 2. For the measured data, the highest mean, median, maximum, standard deviation, and skewness values are recorded at Mykonos, and for the model wave data, the highest values for median, maximum, and standard deviation are encountered at Gulf of Biscay. Gulf of Lion exhibits by far the highest CV P , skewness, and kurtosis values. Clearly, the most energetic wave climate from the examined locations is encountered in the Gulf of Biscay. In Table 3, the values of the main circular statistical parameters (i.e., mean direction m Θ , median direction Θ 0.5 (its interpretation is analogous to the classic median), mean resultant length R Θ (that quantifies the spread of the variable), circular variance V Θ = 1 − R Θ , and circular standard deviation s Θ = −2 ln R Θ , as regards wave direction Θ) are presented. The definitions of the above parameters can be found in [58]. The mean wave direction at Zakynthos and Gulf of Lion is coming from the WSW sector, at Santorini and Gulf of Biscay from NW, at Mykonos from the NNW sector, and at the North Sea location from the Processes 2021, 9, 460 6 of 24 WNW sector. Moreover, the highest and the lowest values of R Θ are encountered in Gulf of Biscay (0.904) and Gulf of Lion (0.130), respectively, while the corresponding values for V Θ can be found in Gulf of Lion (0.870) and Gulf of Biscay (0.096). The value of V Θ for Gulf of Lion, which is close to 1, denotes that the range of the corresponding directions is wide while the value of R Θ for Gulf of Biscay reveals a small dispersion of the data. In Figure 2, the directional histograms of wave energy flux are presented for the examined locations.
sented. The definitions of the above parameters can be found in [58]. The mean wave direction at Zakynthos and Gulf of Lion is coming from the WSW sector, at Santorini and Gulf of Biscay from NW, at Mykonos from the NNW sector, and at the North Sea location from the WNW sector. Moreover, the highest and the lowest values of ̅ are encountered in Gulf of Biscay (0.904) and Gulf of Lion (0.130), respectively, while the corresponding values for can be found in Gulf of Lion (0.870) and Gulf of Biscay (0.096). The value of for Gulf of Lion, which is close to 1, denotes that the range of the corresponding directions is wide while the value of ̅ for Gulf of Biscay reveals a small dispersion of the data. In Figure 2, the directional histograms of wave energy flux are presented for the examined locations.

Parametric Models for Linear Variables
Since there is no theoretical justification for adopting a particular probability distribution for wave energy flux, several distributions are implemented according to specific goodness-of-fit tests. The initial evaluation criteria are based on the performance of two goodness-of-fit tests, the Kolmogorov-Smirnov (K-S) and the Anderson-Darling (A-D) tests that were calculated for 24 conventional distributions for the modelling of ; see also [47,59,60]. The examined distributions are the following: (1) Two-parameter cdfs: Exponential, Gaussian, Rayleigh; The two-and three-parameter pdfs determine the first, second, and third moments of the distribution (i.e., mean value and variance, and skewness), respectively, while the pdfs with four parameters represent additionally kurtosis (fourth moment). Typically, as the number of parameters increases, the distributions tend to be more robust and can provide a better representation of the data. However, distributions with small number of parameters can be more straightforward as regards their definition and use in applications.

Parametric Models for Linear Variables
Since there is no theoretical justification for adopting a particular probability distribution for wave energy flux, several distributions are implemented according to specific goodness-of-fit tests. The initial evaluation criteria are based on the performance of two goodness-of-fit tests, the Kolmogorov-Smirnov (K-S) and the Anderson-Darling (A-D) tests that were calculated for 24 conventional distributions for the modelling of P; see also [47,59,60]. The examined distributions are the following: (1) Two-parameter cdfs: Exponential, Gaussian, Rayleigh; The two-and three-parameter pdfs determine the first, second, and third moments of the distribution (i.e., mean value and variance, and skewness), respectively, while the pdfs with four parameters represent additionally kurtosis (fourth moment). Typically, as the number of parameters increases, the distributions tend to be more robust and can provide a better representation of the data. However, distributions with small number of parameters can be more straightforward as regards their definition and use in applications.
Then, the ten distributions that provided the smallest statistic values of K-S and A-D tests were selected for each examined location. Due to several overlaps between the examined locations and the two tests, the total number of the fitted distributions with a good performance was relatively small; see also Section 6.1. The distributions that proceeded to the next step of evaluation are the following: Burr (BUR), Dagum (DAG), Lognormal (LGN), Fatigue Life (FAL), Pearson 6 (PE6), Log-Logistic (LGL), and Generalized Gamma (GNG). In the Appendix A, the analytic forms of the above-mentioned probability density functions (pdfs) are provided. Moreover, four parametric mixture distributions are also considered: one homogeneous, the Weibull (WW) mixture distribution, and three inhomogeneous, the 2-parameter Weibull-GEV (WGEV), the Exponential-Lognormal (ELGN), and the Exponential-Weibull (EW) mixtures. For the definition of these mixture distributions, see the Appendix A. The parametric mixture distributions were selected due to their suitability in describing right-skewed data and were not considered in the evaluation of the rest conventional distributions by means of K-S and A-D tests. The selected mixture distributions are characterized by rather large values of R 2 a,1 (≥ 0.99). See also Section 5.1.

Parametric Models for Directional Variables
A finite mixture of von Mises (vM) distributions (i.e., the equivalent of the Gaussian distribution for linear variables) has been proved very appropriate for modelling directional variables with more than one mode such as wave or wind direction [40,61]. The pdf of a vM mixture with J-components is defined as follows: where I 0 (κ) is the modified Bessel function of the first kind and zero-th order, i.e., I 0 (κ) = 1 π π 0 e κ cos θ dθ; J is the number of components, κ j and µ j , j = 1, 2, · · · , J, is the concentration parameter and mean direction, respectively, of each individual vM distribution, and ω j 's are the weighting quantities that sum to one.
In order to determine the optimal number of components that should be included in the vM mixture model, the Bayesian information criterion (BIC) is applied. BIC is dependent on the sample size and balances the likelihood and the complexity of the mixture model, i.e., BIC = −2 ln L + p ln n, where L is the maximized value of the likelihood function for the mixture model, p is the number of parameters to be estimated and n is the sample size of each wave dataset.
For various values of j in the range [3, ln n + 1], the optimal number of components for the estimation of the parameters corresponds to the model with the smallest value of BIC. The packages "Directional" [62] and "movMF" [63] were used for the estimation of the BIC values and the implementation of the expectation-maximization algorithm for maximum likelihood estimates of parameters for the mixture model, respectively, available at Comprehensive R Archive Network.

The Johnson-Wehrly Bivariate Model
In this section, the Johnson-Wehrly (JW) bivariate distribution is described for the modelling of linear and directional variables, i.e., wave energy flux and wave direction. The main feature of the JW model is that it is marginally consistent, i.e., the integration of the bivariate model over the support of each of the two variables coincides with the two corresponding marginal pdfs.
Bivariate probability models have been used to describe environmental and metocean parameters, see, e.g., [32,50,[64][65][66][67]. In particular, in [67], the joint pdf of significant wave height H S and mean zero-up crossing wave period T Z was constructed based on the Plackett family for various types of marginal distributions while [32] compared fits

Univariate Distributions
As mentioned in Section 3.1, the preliminary selection of the most efficient distributions for modelling wave energy flux was based on the K-S and A-D tests. Since the two tests are not equivalent, with the latter being more sensitive to deviations in the tails of the distribution, they result, in general, in different "best-fit" distributions; for this reason, the first ten best-fit models provided by each test are selected. However, the total number of the pre-final distributions is rather small due to many overlaps between the two tests.
The definitive evaluation of the univariate distributions is based on the coefficient of determination R 2 a,1 , adopted in similar studies such as wind speed modelling; see, e.g., [47,[69][70][71][72]. R 2 a,1 measures the degree of association between the estimated (predicted) and observed cumulative probabilities of a selected distribution and is defined as follows: whereF i is the estimated cumulative distribution (obtained from each distribution examined), F i is the empirical (observed) distribution function and F = (1/n) ∑ n i=1 F i . A value of R 2 a,1 close to 1 indicates that the candidate distribution provides a good fit to the observed (measured) data.

Bivariate Distributions
The evaluation of the bivariate distribution fit is based on five bin-specific measures of performance, namely, the root mean square error (RMSE), mean absolute error (MAE), adjusted coefficient of determination R 2 a,2 , index of agreement (IA), and χ 2 statistic. See Table 4 and [47] for the definitions of the above goodness-of-fit measures. These global measures are derived from various combinations (scales) of the fundamental quantity ij denote the observed and the estimated fraction of points that belong to the (i, j)-cell, respectively. Evidently, ∆ o,e ij indicates the magnitude of the difference between the observed and the estimated frequency of occurrence for a specific bin. Nevertheless, bin-specific measures suffer from the subjective selection of the bin size; see also Section 6.5 in [47], for a discussion on this issue. In the present work, the selected bin size for wave energy flux is 1 kW/m and for wave direction is 5 deg. A "good model" is characterized by small values (as close to 0 as possible) of MAE, RMSE, and χ 2 , and high values (as close to 1 as possible) of IA and R 2 a,2 .
ij ): difference between the observed (estimated) and the mean frequency, respectively; n: sample size; N B : total number of non-empty bins; N B : total number of bins; q: number of parameters estimated for the particular distribution.
Since ambiguity can be observed between the results of the above measures of performance (i.e., when the values of the above goodness-of-fit measures do not converge to the same model), a final, very simple, synthetic criterion (SC) is implemented that combines all the above measures into a single expression. According to [49], this criterion can be defined as follows: SC is defined as the ratio of various goodness-of-fit statistics, where the nominator corresponds to statistics that suggest a good fit when they take high values (in our case R 2 a,2 and I A), and the denominator corresponds to statistics that suggest a good fit when they take low values (in our case RMSE, χ 2 , and MAE). Therefore, in this case, large values of SC suggest a good fit. The nominator of Equation (8) is in the range [0, 1] while the denominator in [ 0, ∞) . Let it also be noted that since there is not an established acceptable threshold value for any of the abovementioned criteria (including SC) clearly dictating that a distribution, say F, is a "good" distribution for modelling the data, they should be interpreted and applied as comparative tools and not as absolute decisionmaking measures.
Other approaches for combining different goodness-of-fit measures in order to select the best-fit distribution are presented in [73], where the data envelopment analysis is introduced, and [74], where a weighted sum of the evaluation criteria is implemented.  In summary, the following univariate conventional distributions are further assessed in the bivariate case (as marginal cdfs for P): The parameters for each candidate conventional distribution, along with the corresponding parameters and the weighting parameter ω of the two mixture distributions, are presented in Table 5. The presented parameters are the following:
In Figure 3, the best-fit distributions along with the corresponding histograms for wave power flux at the examined locations are presented. The range of the values on the x-axis varies, depending on the corresponding behaviour at each examined location, while the bin size of the histograms was kept the same, i.e., 1 kW/m.

Univariate Distributions for Mean Wave Direction
In Table 6, the parameters of the mixtures of vM distributions for Θ modelling are summarized. In terms of BIC, wave direction at Athos, Zakynthos, and North Sea is best described with five components of the vM mixtures, at Mykonos with four, at Santorini with seven, and at the Gulf of Lion with twelve.
In Figure 4, the histograms of mean wave direction for the examined locations are shown along with the plotted lines for vM mixtures. −axis varies, depending on the corresponding behaviour at each examined location, while the bin size of the histograms was kept the same, i.e., 1 kW/m.

Univariate Distributions for Mean Wave Direction
In Table 6, the parameters of the mixtures of vM distributions for modelling are summarized. In terms of BIC, wave direction at Athos, Zakynthos, and North Sea is best described with five components of the vM mixtures, at Mykonos with four, at Santorini with seven, and at the Gulf of Lion with twelve.   In Figure 4, the histograms of mean wave direction for the examined locations are shown along with the plotted lines for vM mixtures.

Joint Pdf for Wave Energy Flux and Direction
The parameters of the vM mixture ( , , ) for modelling are also estimated (not provided here) for the examined locations and pdfs of wave energy flux. Since is a smooth function, two components of the vM mixture are adequate for its description.

Joint Pdf for Wave Energy Flux and Direction
The parameters of the vM mixture (κ, µ, ω) for modelling f Ψ PΘ are also estimated (not provided here) for the examined locations and pdfs of wave energy flux. Since f Ψ PΘ is a smooth function, two components of the vM mixture are adequate for its description.
In Table 7, the values of the goodness-of-fit criteria for the bivariate parametric distributions are provided for the examined locations and datasets. Boldface numbers denote the best value obtained from all the examined bivariate families. For Mykonos, the JW model with the GNG marginal distribution for wave flux, provides the optimum values for the FC criterion. Let it be noted that for Mykonos, the optimal marginal distribution according to R 2 a,1 , was WGEV mixture, while the GNG distribution provided the third best fit. For Santorini and Zakynthos, the ING distribution provides the optimum values for the FC criterion, although the optimal marginal distribution was LGN distribution. It is reminded that the ING distribution provided the second-best fit for Santorini and the third best fit for Zakynthos in the univariate case.
For Gulf of Lion and North Sea, the JW model with the ELGN marginal distribution (for wave flux) provides the optimum FC value. The same distribution provided the best fit for the univariate case for both locations. Finally, for Gulf of Biscay, the JW model with the WGEV mixture marginal distribution provides the optimum values for the FC criterion. Let it be noted that the optimal marginal distribution according to R 2 a,1 was LGN distribution while the WGEV mixture distribution provided the third best fit for wave energy flux.
Regarding the involvement of the univariate distributions per se in the bivariate models, mixture distributions (ELGN and WGEV) for wave energy flux provide the best bivariate fits for the locations where numerical wave model data are available. For Santorini and Zakynthos, ING distribution provides the best bivariate fits while for Mykonos GNG provide the best bivariate fit. It can also be deduced that the best univariate conventional distribution (in terms of R 2 a,1 ) for wave energy flux does not necessarily provides the best fit in the bivariate case.
In Figure 5, the best bivariate model for modelling the variables of interest is plotted for each location.

Application-Oriented Results
The preliminary assessment of the available wave energy in an area requires the estimation of the wave energy flux per sea state that is a function of significant wave height and energy period (see Equation (3)). Due to the temporal and spatial variability of , , wave direction , and sea water density , more detailed information as regards wave flux for different occurrences of these variables at the specific area is required during the stage of wave farm design and layout. In this respect, for the evaluation of different aspects of the wave energy flux at a particular location, the wave energy flux-wave direction distribution should be used.
For a series of sea states , , , , , = 1,2, ⋯ , , in an area, and for constant seawater density, the corresponding realizations , of the bivariate random variable ( , ) can be easily obtained, and its corresponding pdf , ( , ) can be estimated following the procedure described in Section 4. The mean wave energy flux for all wave directions can be obtained as follows:

Application-Oriented Results
The preliminary assessment of the available wave energy in an area requires the estimation of the wave energy flux per sea state that is a function of significant wave height H S and energy period T e (see Equation (3)). Due to the temporal and spatial variability of H S , T e , wave direction Θ, and sea water density ρ, more detailed information as regards wave flux for different occurrences of these variables at the specific area is required during the stage of wave farm design and layout. In this respect, for the evaluation of different aspects of the wave energy flux at a particular location, the wave energy flux-wave direction distribution should be used.
For a series of sea states H S,i , T e,i , Θ i , i = 1, 2, · · · , n, in an area, and for constant seawater density, the corresponding realizations P i , Θ i of the bivariate random variable (P, Θ) can be easily obtained, and its corresponding pdf f P,Θ (p, θ) can be estimated following the procedure described in Section 4. The mean wave energy flux for all wave directions can be obtained as follows: where P min , P max are the minimum and maximum values of P observed in the area, respectively (or, equivalently, of 0.49H 2 S T e ). From Equation (9), useful and detailed characteristics of the available directional wave energy flux in an area can be obtained. These characteristics can be further used as design parameters for the WEC and the wave farm layout.
For example, the conditional mean value of P with respect to Θ can be derived as follows: where f P|Θ (p|θ ) is the conditional pdf that is defined as follows: From the above relations, the wave energy flux per wave direction for various ranges of values of H 2 S T e , the wave energy flux per 0.49H 2 S T e for various values of Θ, and the distribution of wave energy flux for some critical directions of interest can also be obtained.
In Figure 6, the expected values of wave energy flux are depicted for the entire range of wave directions in the examined areas. As shown in Figure where , are the minimum and maximum values of observed in the area, respectively (or, equivalently, of 0.49 2 ). From Equation (9), useful and detailed characteristics of the available directional wave energy flux in an area can be obtained. These characteristics can be further used as design parameters for the WEC and the wave farm layout.
For example, the conditional mean value of with respect to can be derived as follows: where | ( | ) is the conditional pdf that is defined as follows: From the above relations, the wave energy flux per wave direction for various ranges of values of 2 , the wave energy flux per 0.49 2 for various values of , and the distribution of wave energy flux for some critical directions of interest can also be obtained.
In Figure 6, the expected values of wave energy flux are depicted for the entire range of wave directions in the examined areas. As shown in Figure 6, for Zakynthos, Mykonos, and Gulf of Biscay, the highest values of the expected values of wave energy flux are concentrated in relatively narrow bands of wave direction. For Santorini, North Sea, and Gulf of Lion, the wave direction sectors that provide large values for the expected wave energy flux are rather wide. In Figure 7, the distribution of wave energy flux for particular directions of interest, i.e., the mean and the most probable direction (based on the histograms of Figure 4) are depicted for the examined locations. It is observed that wave energy flux has the highest likelihood for the most probable wave direction at all examined locations. Let it be noted that regarding the examined locations with measured data, Mykonos is the only one where the likelihood remains high enough for values of wave energy flux above 5 kW/m, denoting that this area can be candidate for a more detailed assessment of wave energy resource. Moreover, it is clear that the Gulf of Biscay is characterized by a very energetic wave climate (especially for the most probable wave direction), which renders this location very promising for wave energy utilization. In Figure 7, the distribution of wave energy flux for particular directions of interest, i.e., the mean and the most probable direction (based on the histograms of Figure 4) are depicted for the examined locations. It is observed that wave energy flux has the highest likelihood for the most probable wave direction at all examined locations. Let it be noted that regarding the examined locations with measured data, Mykonos is the only one where the likelihood remains high enough for values of wave energy flux above 5 kW/m, denoting that this area can be candidate for a more detailed assessment of wave energy resource. Moreover, it is clear that the Gulf of Biscay is characterized by a very energetic wave climate (especially for the most probable wave direction), which renders this location very promising for wave energy utilization.

Conclusions
The joint modelling of wave energy flux and wave direction in an area is important for wave energy projects. For instance, wave resource assessment issues and, in turn, WEC siting problems can be mitigated by constructing the most appropriate bivariate model based on real world data, attaining in this way economic viability of the corresponding projects. In this work, the Johnson-Wehrly parametric model was examined for the joint probabilistic description of wave energy flux and wave direction, after evaluating various conventional distributions for the linear variable. The proposed parametric bivariate models were constructed by considering the corresponding marginal distributions along with an appropriate dependence structure of the involved variables. The performance of the bivariate models was evaluated based on a global goodness-of-fit criterion, formed as a combination of RMSE, MAE, IA, 2 , and ,2 2 . The models were applied for six long-term wave energy flux and wave direction time series, which were obtained from in situ measurements in the Greek Seas and wave model results in the western Mediterranean Sea, Gulf of Biscay, and the North Sea.
From the results, it was found that as regards wave energy flux, wave measurements were best modelled by the Inverse Gaussian and Generalized Gamma distributions while Exponential-Lognormal and Weibull-GEV mixture distributions provided overall the best fits for the examined locations where numerical model results were available. However, the above distributions were not the optimum ones in the evaluation of the univariate case, highlighting thus the necessity to evaluate a wide range of models for wave energy applications. The best bivariate model was also applied for the estimation of wave energy flux per wave direction as a real-world scenario in wave farm design for the examined locations.
This research revealed that an important prerequisite related to WEC developers is the provision of the corresponding power matrices for different mean wave directions. In this way, uncertainties as regards the annual wave energy production would be reduced while the optimization and the siting of a wave farm can be greatly facilitated.

Conclusions
The joint modelling of wave energy flux and wave direction in an area is important for wave energy projects. For instance, wave resource assessment issues and, in turn, WEC siting problems can be mitigated by constructing the most appropriate bivariate model based on real world data, attaining in this way economic viability of the corresponding projects. In this work, the Johnson-Wehrly parametric model was examined for the joint probabilistic description of wave energy flux and wave direction, after evaluating various conventional distributions for the linear variable. The proposed parametric bivariate models were constructed by considering the corresponding marginal distributions along with an appropriate dependence structure of the involved variables. The performance of the bivariate models was evaluated based on a global goodness-of-fit criterion, formed as a combination of RMSE, MAE, IA, χ 2 , and R 2 a,2 . The models were applied for six longterm wave energy flux and wave direction time series, which were obtained from in situ measurements in the Greek Seas and wave model results in the western Mediterranean Sea, Gulf of Biscay, and the North Sea.
From the results, it was found that as regards wave energy flux, wave measurements were best modelled by the Inverse Gaussian and Generalized Gamma distributions while Exponential-Lognormal and Weibull-GEV mixture distributions provided overall the best fits for the examined locations where numerical model results were available. However, the above distributions were not the optimum ones in the evaluation of the univariate case, highlighting thus the necessity to evaluate a wide range of models for wave energy applications. The best bivariate model was also applied for the estimation of wave energy flux per wave direction as a real-world scenario in wave farm design for the examined locations.
This research revealed that an important prerequisite related to WEC developers is the provision of the corresponding power matrices for different mean wave directions. In this way, uncertainties as regards the annual wave energy production would be reduced while the optimization and the siting of a wave farm can be greatly facilitated. Acknowledgments: Flora Karathanasi was partially supported for this work, which is part of a project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 857586.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
English letters F X (x) cumulative distribution function of the random variable X f X (x) probability density function of the random variable X F XY (x, y) joint (bivariate) cumulative distribution function of the random variables X, Y f XY (x, y) joint (bivariate) probability density function of the random variables X, Y F X estimate of F X obtained from a fitted analytic model I 0 modified Bessel function of the first kind and zero order n sample size N B total number of bins N B total number of non-empty bins f frequency/probability of occurrence R 2 a,1 coefficient of determination for the univariate case R 2 a,2 coefficient of determination for the bivariate case X linear random variable x(F X ) quantile function of the random variable X (inverse cumulative distribution function) P W , P wave energy flux H S significant wave height T e mean wave energy period T Z mean zero-upcrossing wave period r PΘ linear-circular correlation coefficient between the linear variable P and the circular variable Θ Greek letters Θ W , Θ mean wave direction (random variable) Φ(·), ϕ(·) Gaussian cumulative distribution and probability density functions, respectively ψ PΘ "correlation-type" parameters for linear P and directional Θ random variables in the Johnson-Wehrly model where Φ(·) is the Laplace integral: x 0 e −t 2 /2 dt.