Idealized models of the joint probability distribution of wind speeds

The joint probability distribution of wind speeds at two separate locations in space or points in time completely characterizes the statistical dependence of these two quantities, providing more information than linear measures such as correlation. In this study, we consider two models of the joint distribution of wind speeds obtained from idealized models of the dependence structure of the horizontal wind velocity components. The bivariate Rice distribution follows from assuming that the wind components have Gaussian and isotropic fluctuations. The bivariate Weibull distribution arises from power law transformations of wind speeds corresponding to vector components with Gaussian, isotropic, mean-zero variability. Maximum likelihood estimates of these distributions are compared using wind speed data from the mid-troposphere, from different altitudes at the Cabauw tower in the Netherlands, and from scatterometer observations over the sea surface. While the bivariate Rice distribution is more flexible and can represent a broader class of dependence structures, the bivariate Weibull distribution is mathematically simpler and may be more convenient in many applications. The complexity of the mathematical expressions obtained for the joint distributions suggests that the development of explicit functional forms for multivariate speed distributions from distributions of the components will not be practical for more complicated dependence structure or more than two speed variables.

1 Introduction

A fundamental issue in the characterization of atmospheric variability is that of dependence: how the state of one atmospheric variable is related to that of another at a different position in space, or point in time. The simplest measure of statistical dependence, the correlation coefficient, is a natural measure for Gaussian-distributed quantities but does not fully characterize dependence for non-Gaussian variables. The most general representation of dependence between two or more quantities is their joint probability distribution. The joint probability distribution for a multivariate Gaussian is well known, and expressed in terms of the mean and covariance matrix (e.g. Wilks, 2005; von Storch and Zwiers, 1999). No such general expressions for non-Gaussian multivariate distributions exist. Copula theory (e.g. Schlözel and Friederichs, 2008) allows joint distributions to be constructed from specified marginal distributions. However, which copula model to use for a given analysis is not generally known a priori and is usually determined empirically through a statistical model selection exercise.

The present study considers the bivariate joint probability distribution of wind speeds. As these are quantities which are by definition bounded below by zero, the joint distribution and the marginal distributions are non-Gaussian. While the correlation structure (equivalently, the power spectrum) of wind speeds in time (e.g. Brown et al., 1984; Schlax et al., 2001; Gille, 2005; Monahan, 2012b) and in space (e.g. Carlin and Haslett, 1982; Nastrom and Gage, 1985; Wylie et al., 1985; Brown and Swail, 1988; Xu et al., 2011) has been considered, relatively little attention has been paid to developing expressions for the joint distribution. Previous studies have used copula methods to model horizontal spatial dependence of wind speeds for wind power applications (Grothe and Schnieders, 2011; Louie, 2012; Veeramachaneni et al., 2015) and dependence of daily wind speed maxima (Schlözel and Friederichs, 2008). While these earlier analyses have focused on probabilistic modelling of simultaneous wind speed values at different spatial locations in the horizontal, dependence structures in the vertical (e.g. for vertical interpolation of wind speeds) or in time are also of interest. For example, an analysis in which the need for an explicit parametric form for the joint distribution of wind speeds at different altitudes has arisen is the hidden Markov model (HMM) analysis considered in Monahan et al. (2015). In an HMM analysis of continuous variables, it is necessary to specify the parametric form of the joint distribution within each hidden state. In Monahan et al. (2015), the joint distribution of wind speeds at 10 and 200 m and the potential temperature difference between these altitudes was modelled as a multivariate Gaussian distribution, despite the fact that for at least the speeds this distribution cannot be strictly correct. This pragmatic modelling approximation was made because of the absence of a more appropriate parametric distribution for the quantities being considered. The alternative approach of using the wind components at the two levels (for which the multivariate Gaussian model may be a better approximation) instead of the speeds directly has the downside of increasing the dimensionality of the state vector from three to five, dramatically increasing the number of parameters to be estimated (with the covariance matrices in particular increasing from 9 to 25 elements) and reducing the statistical robustness of the results.

A number of previous studies have constructed univariate speed distributions starting from models for the joint distribution of the horizontal components (e.g. Cakmur et al., 2004; Monahan, 2007; Drobinski et al., 2015). One useful benefit of this approach is that it allows the statistics of the speed and the components to be related to each other. The specific goal of the present study is to extend this approach to bivariate distributions, constructing models of the joint probability distribution of wind speeds that are directly connected to the joint distributions of the horizontal components of the wind. As the following results will demonstrate, generalizing this approach to the bivariate speed distribution results in rather complicated mathematical expressions. Expressions for multivariate distributions of more than two speeds will be even more complicated, and may not be analytically tractable. Through the analysis of the bivariate speed distribution, we will probe how far the development of closed-form, analytic expressions for parametric speed distributions based on distributions of components can practically be extended.

Both Weibull and Rice distributions have been used to model the univariate wind speed distribution (cf. Carta et al., 2009; Monahan, 2014; Drobinski et al., 2015), and models of multivariate distributions with Weibull or Ricean marginals have been developed (e.g. Crowder, 1989; Lu and Bhattacharyya, 1990; Kotz et al., 2000; Sagias and Karagiannidis, 2005; Yacoub et al., 2005; Mendes and Yacoub, 2007; Villanueva et al., 2013). Much of this work has been done in the context of wireless communications (Sagias and Karagiannidis, 2005; Yacoub et al., 2005; Mendes and Yacoub, 2007): the present study builds upon the results of these earlier analyses.

The two probability density functions (pdfs) we will consider, the bivariate Rice and Weibull distributions, both start with simple assumptions regarding the distributions of the wind components. The bivariate Rice distribution follows directly from the assumption of Gaussian components with isotropic variance, but nonzero mean. In contrast, the bivariate Weibull distribution is obtained from nonlinear transformations of the magnitudes of Gaussian, isotropic, mean-zero components. While the univariate Weibull distribution has been found to generally be a better fit to observed wind speeds than the univariate Rice distribution (particularly over the oceans, e.g. Monahan, 2006, 2007), the direct connection of the Rice distribution to the distribution of the components (which the Weibull distribution does not have) is useful from a modelling and theoretical perspective (e.g. Cakmur et al., 2004; Monahan, 2012a; Culver and Monahan, 2013; Sun and Monahan, 2013; Drobinski et al., 2015). The six-parameter bivariate Rice distribution that we will consider is more flexible than the five-parameter bivariate Weibull distribution, and able to model a broader range of dependence structures. Furthermore, it is directly connected to the univariate distributions and dependence structure of the wind components. However, the bivariate Weibull distribution is mathematically much simpler than the bivariate Rice distribution and easier to use in practice. Other flexible bivariate distributions for non-negative random variables exist, such as the α-μ distribution discussed in Yacoub (2007). Because the Weibull and Rice distributions are common models for the univariate wind speed distributions, this study will focus specifically on their bivariate generalizations.

The bivariate Rice and Weibull distributions are developed in Sect. 2, starting from discussion of the bivariate Rayleigh distribution (which is a limiting case of both of the other models). In this section, we repeat some of the formulae obtained by Sagias and Karagiannidis (2005), Yacoub et al. (2005), and Mendes and Yacoub (2007) for completeness and because of notational differences between this study and the earlier ones. In Sect. 3, the ability of these distributions to model wind speed data from the middle troposphere, and from the near-surface flow over land and the ocean, is considered. Examples of dependence structures in both space (horizontally and vertically) and in time are considered. A discussion and conclusions are presented in Sect. 4.

2 Empirical models of the bivariate wind speed distribution

As a starting point for developing models of the bivariate wind speed distribution, we consider the joint distribution of the horizontal wind vector components ${\mathit{u}}_i=(u_i,v_i)$, i=1, 2 (where the subscripts i denote wind vectors at two different locations, two different points in time, or both). In particular, we assume that

1. the two orthogonal wind components are marginally Gaussian with isotropic and uncorrelated fluctuations:

   $$\begin{array}{}\text{(1)}& \left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)\sim \mathcal{N}\left[\left(\begin{array}{c}{\stackrel{\mathrm{‾}}{u}}_i\\ {\stackrel{\mathrm{‾}}{v}}_i\end{array}\right),\left(\begin{array}{cc}{\mathit{\sigma }}_i^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& {\mathit{\sigma }}_i^{\mathrm{2}}\end{array}\right)\right]i=\mathrm{1},\mathrm{2},\end{array}$$

2. and the cross-correlation matrix of the two vectors is

   $$\begin{array}{ll}\text{(2)}& {\displaystyle }\mathrm{corr}({\mathit{u}}_{\mathrm{1}},{\mathit{u}}_{\mathrm{2}})& {\displaystyle }=\left(\begin{array}{cc}\mathrm{corr}(u_{\mathrm{1}},u_{\mathrm{2}})& \mathrm{corr}(u_{\mathrm{1}},v_{\mathrm{2}})\\ \mathrm{corr}(v_{\mathrm{1}},u_{\mathrm{2}})& \mathrm{corr}(v_{\mathrm{1}},v_{\mathrm{2}})\end{array}\right){\displaystyle }& {\displaystyle }=\left(\begin{array}{cc}{\mathit{\mu }}_{\mathrm{1}}& {\mathit{\mu }}_{\mathrm{2}}\\ -{\mathit{\mu }}_{\mathrm{2}}& {\mathit{\mu }}_{\mathrm{1}}\end{array}\right)=\left(\begin{array}{cc}\mathit{\rho }\cos \mathit{\psi }& \mathit{\rho }\sin \mathit{\psi }\\ -\mathit{\rho }\sin \mathit{\psi }& \mathit{\rho }\cos \mathit{\psi }\end{array}\right),\end{array}$$

   where we have expressed the correlations in both Cartesian and polar coordinates: $({\mathit{\mu }}_{\mathrm{1}},{\mathit{\mu }}_{\mathrm{2}})=(\mathit{\rho }\cos \mathit{\psi },\mathit{\rho }\sin \mathit{\psi })$, with $\mathrm{0}\le \mathit{\rho }={\left({\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\le \mathrm{1}$. This assumed correlation structure implies that the correlation matrix becomes diagonal when the vector u₂ is rotated through the angle −ψ:

   $$\begin{array}{}\text{(3)}& \mathrm{corr}({\mathit{u}}_{\mathrm{1}},R(-\mathit{\psi }\left){\mathit{u}}_{\mathrm{2}}\right)=\left(\begin{array}{cc}\mathit{\rho }& \mathrm{0}\\ \mathrm{0}& \mathit{\rho }\end{array}\right),\end{array}$$

   where

   $$\begin{array}{ll}{\displaystyle }R(-\mathit{\psi }){\mathit{u}}_{\mathrm{2}}& {\displaystyle }=\left(\begin{array}{cc}\cos \mathit{\psi }& \sin \mathit{\psi }\\ -\sin \mathit{\psi }& \cos \mathit{\psi }\end{array}\right)\left(\begin{array}{c}{u}_{\mathrm{2}}\\ v_{\mathrm{2}}\end{array}\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }=\left(\begin{array}{c}{u}_{\mathrm{2}}\cos \mathit{\psi }+v_{\mathrm{2}}\sin \mathit{\psi }\\ -u_{\mathrm{2}}\sin \mathit{\psi }+v_{\mathrm{2}}\cos \mathit{\psi }\end{array}\right).\end{array}$$

The joint distribution of the horizontal components resulting from these assumptions is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(u_{\mathrm{1}},u_{\mathrm{2}},v_{\mathrm{1}},v_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{(u_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{(v_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right.\right.\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(u_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}+{\displaystyle \frac{(v_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}\left[\right(u_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}\left)\right(u_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}})+(v_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}\left)\right(v_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}\left)\right]}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }\left.\left.-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}\left[\left(u_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}\right)(v_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}})-(v_{\mathrm{1}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}})(u_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}})\right]}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\right)\end{array}$$

(Mendes and Yacoub, 2007).

Note that only considering the horizontal components of the wind vector implicitly restricts the resulting distributions to timescales sufficiently long that the vertical component of the wind contributes negligibly to the speed.

2.1 Bivariate Rayleigh distribution

The joint distributions of the speeds $w_i=\sqrt{u_i^{\mathrm{2}}+v_i^{\mathrm{2}}}$ obtained from the pdf Eq. (5) when both vector wind components are mean-zero is the bivariate Rayleigh distribution (e.g. Battjes, 1969). Transforming variables to wind speed w_i and direction ${\mathit{\theta }}_i={\tan }^{-\mathrm{1}}(v_i/u_i)$, the joint distribution (Eq. 5) with ${\stackrel{\mathrm{‾}}{u}}_i=\stackrel{\mathrm{‾}}{v_i}=\mathrm{0}$, $i=\mathrm{1},\mathrm{2}$, becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}},{\mathit{\theta }}_{\mathrm{1}},{\mathit{\theta }}_{\mathrm{2}})=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{\left(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{w_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right]\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }\times \exp \left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\cos ({\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}+\mathit{\psi })\right).\end{array}$$

Integrating over the wind directions to obtain the marginal distribution for the wind speeds, we obtain

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}})=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{w_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right]\right)\\ \text{(7)}& {\displaystyle }& {\displaystyle }{I}_o\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right),\end{array}$$

where we have used the fact that

$$\begin{array}{}\text{(8)}& \underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{e}^{\mathit{\alpha }\cos \mathit{\theta }}\cos k\mathit{\theta }\mathrm{d}\mathit{\theta }=\mathrm{2}\mathit{\pi }{I}_k\left(\mathit{\alpha }\right),\end{array}$$

where I_k(z) is the modified Bessel function of order k. Note that the correlation angle ψ drops out after integration over θ₁ and θ₂. As a result, for the bivariate Rayleigh distribution, p(w₁,w₂) depends only on the three parameters $({\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}},\mathit{\rho })$.

For ρ=0, p(w₁,w₂) factors as the product of the marginal distributions of w₁ and w₂:

$$\begin{array}{}\text{(9)}& p(w_{\mathrm{1}},w_{\mathrm{2}})=\left[{\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)\right]\left[{\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)\right],\end{array}$$

and the wind speeds are independent. As ρ→1, we can use the asymptotic result

$$\begin{array}{}\text{(10)}& I_{\mathrm{0}}\left(x\right)\sim {\displaystyle \frac{e^x}{\sqrt{\mathrm{2}\mathit{\pi }x}}}(x\gg \mathrm{1})\end{array}$$

to show that

$$\begin{array}{}\text{(11)}& p(w_{\mathrm{1}},w_{\mathrm{2}})\to {\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)\mathit{\delta }\left({\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}-{\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\right),\end{array}$$

where δ(⋅) is the Dirac delta function. In this limit, w₁ and w₂ are perfectly correlated and Rayleigh distributed.

Figure 1 Example bivariate Rayleigh distributions p(w₁,w₂) for $\mathit{\rho }=\mathrm{0},\mathrm{0.5}$, and 0.85, with w₁ and w₂ scaled respectively by σ₁ and σ₂. The upper and left subpanels show the marginal distributions of w₁ and w₂ respectively. These marginal distributions are the same for all three panels.

Moments of the bivariate Rayleigh distribution are given by

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }\mathsf{E}\left\{w_{\mathrm{1}}^m{w}_{\mathrm{2}}^n\right\}={\displaystyle }& {\displaystyle }{\mathrm{2}}^{m/\mathrm{2}}{\mathrm{2}}^{n/\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}^m{\mathit{\sigma }}_{\mathrm{2}}^n\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{m}{\mathrm{2}}}\right)\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{n}{\mathrm{2}}}\right){}_{\mathrm{2}}{F}_{\mathrm{1}}\left(-{\displaystyle \frac{m}{\mathrm{2}}},-{\displaystyle \frac{n}{\mathrm{2}}},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}}\right),\end{array}$$

where ${}_{\mathrm{2}}{F}_{\mathrm{1}}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma };z)$ is the hypergeometric function (Gradshteyn and Ryzhik, 2000). In particular, we have

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{mean}\left(w_i\right)={\mathit{\sigma }}_i\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}},\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{var}\left(w_i\right)=\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right),\text{(15)}& {\displaystyle }& {\displaystyle }\mathrm{corr}(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}-\mathit{\pi }}}\left({}_{\mathrm{2}}{F}_{\mathrm{1}}\left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}}\right]-\mathrm{1}\right).\end{array}$$

Because ${}_{\mathrm{2}}{F}_{\mathrm{1}}(-\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}})$ is an increasing function of ρ² with ${}_{\mathrm{2}}{F}_{\mathrm{1}}(-\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2},\mathrm{1};\mathrm{0})=\mathrm{1}$, corr(w₁,w₂) must be non-negative for the bivariate Rayleigh distribution.

Plots of p(w₁,w₂) for the three values ρ=0, 0.5, and 0.85 are presented in Fig. 1, along with the marginal distributions of w₁ and w₂ (which are the same for all three panels). The marginal distributions are positively skewed and the contours of the joint distributions are more tightly concentrated below and to the left of their peaks than elsewhere. As expected, the distributions become more tightly concentrated around the 1 : 1 line as the dependence parameter ρ increases.

2.2 Bivariate Rice distribution

The assumptions leading to the bivariate Rayleigh distribution are too restrictive to model observed wind speeds in most circumstances. A more general distribution results from assuming that the wind components are Gaussian, isotropic, and uncorrelated, but with nonzero mean (Eq. 5).

Changing variables to wind speed w_i and direction θ_i, the joint distribution becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}},{\mathit{\theta }}_{\mathrm{1}},{\mathit{\theta }}_{\mathrm{2}})={\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{\left(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\right(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{(w_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}^{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{(w_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}^{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right.\right.\\ {\displaystyle }& {\displaystyle }\left.\left.-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}({\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}({\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\right)\\ {\displaystyle }& {\displaystyle }\times \exp \left(\right[{\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}{a}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}+{\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}{a}_{\mathrm{2}}\cos {\mathit{\theta }}_{\mathrm{2}}+{\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}{b}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}+{\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}{b}_{\mathrm{2}}\sin {\mathit{\theta }}_{\mathrm{2}}\\ \text{(16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\left({\mathit{\mu }}_{\mathrm{1}}\cos ({\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})+{\mathit{\mu }}_{\mathrm{2}}\sin ({\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\right)\left]\right),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }{a}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}-{\mathit{\mu }}_{\mathrm{1}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}-{\mathit{\mu }}_{\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\right)\\ \text{(17)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\cos {\mathit{\varphi }}_{\mathrm{1}}-\mathit{\rho }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\cos ({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\psi })\right],{\displaystyle }{b}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}-{\mathit{\mu }}_{\mathrm{1}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}+{\mathit{\mu }}_{\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathcal{U}}_{\mathrm{1}}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\sin {\mathit{\varphi }}_{\mathrm{1}}-\mathit{\rho }{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathcal{U}}_{\mathrm{2}}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\sin ({\mathit{\varphi }}_{\mathrm{2}}-\mathit{\psi })\right],{\displaystyle }{a}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}-{\mathit{\mu }}_{\mathrm{1}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}+{\mathit{\mu }}_{\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\right)\\ \text{(19)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\cos {\mathit{\varphi }}_{\mathrm{2}}-\mathit{\rho }{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\cos ({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\psi })\right],{\displaystyle }{b}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}-{\mathit{\mu }}_{\mathrm{1}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}-{\mathit{\mu }}_{\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\right)\\ \text{(20)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathcal{U}}_{\mathrm{2}}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\sin {\mathit{\varphi }}_{\mathrm{2}}-\mathit{\rho }{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathcal{U}}_{\mathrm{1}}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\sin ({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\psi })\right],\end{array}$$

and where we have defined the magnitude and direction of the mean vector wind:

$$\begin{array}{}\text{(21)}& ({\stackrel{\mathrm{‾}}{u}}_i,{\stackrel{\mathrm{‾}}{v}}_i)={\stackrel{\mathrm{‾}}{\mathcal{U}}}_i(\cos {\mathit{\varphi }}_i,\sin {\mathit{\varphi }}_i).\end{array}$$

The marginal distribution for the wind speeds is obtained by integrating the joint distribution over θ₁ and θ₂. To evaluate this integral, we make use of the result

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\exp \left[{\mathit{\alpha }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}+{\mathit{\alpha }}_{\mathrm{2}}\cos {\mathit{\theta }}_{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}}\sin {\mathit{\theta }}_{\mathrm{2}}\right.\\ {\displaystyle }& {\displaystyle }\left.+\mathit{\gamma }\cos ({\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}+\mathit{\psi })\right]\mathrm{d}{\mathit{\theta }}_{\mathrm{1}}\mathrm{d}{\mathit{\theta }}_{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{ϵ}}_k\cos \left[k\left({\tan }^{-\mathrm{1}}{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{1}}}}-{\tan }^{-\mathrm{1}}{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2}}}{{\mathit{\alpha }}_{\mathrm{2}}}}+\mathit{\psi }\right)\right]\\ \text{(22)}& {\displaystyle }& {\displaystyle }\times I_k\left(\sqrt{{\mathit{\alpha }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{1}}^{\mathrm{2}}}\right)I_k\left(\sqrt{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{2}}^{\mathrm{2}}}\right)I_k\left(\mathit{\gamma }\right),\end{array}$$

where

$$\begin{array}{}\text{(23)}& {\mathit{ϵ}}_k=\left\{\begin{array}{cc}\mathrm{1}& k=\mathrm{0}\\ \mathrm{2}& k\ne \mathrm{0}\end{array}\right.\end{array}$$

and it is important that ${\tan }^{-\mathrm{1}}(b/a)$ be evaluated as the angle between the vector (a,b) and the vector (1,0) (that is, as the four-quadrant inverse tangent). Equation (22) follows from the Fourier series

$$\begin{array}{}\text{(24)}& e^{c\cos \mathit{\theta }}=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{ϵ}}_k{I}_k\left(c\right)\cos k\mathit{\theta }\end{array}$$

along with repeated use of trigonometric identities and the integral Eq. (8).

Finally, we obtain the expression for the bivariate Rice distribution (Mendes and Yacoub, 2007)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{(w_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}^{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{(w_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}^{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right.\right.\\ {\displaystyle }& {\displaystyle }\left.\left.-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}({\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}({\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}-{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}})}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\right)\\ {\displaystyle }& {\displaystyle }\times \sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{ϵ}}_k\cos \left[k\left({\tan }^{-\mathrm{1}}{\displaystyle \frac{b_{\mathrm{1}}}{a_{\mathrm{1}}}}-{\tan }^{-\mathrm{1}}{\displaystyle \frac{b_{\mathrm{2}}}{a_{\mathrm{2}}}}+{\tan }^{-\mathrm{1}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{1}}}}\right)\right]\\ {\displaystyle }& {\displaystyle }\times I_k\left({\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\sqrt{a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}}}\right)I_k\left({\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\sqrt{a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}}}\right)\\ \text{(25)}& {\displaystyle }& {\displaystyle }\times I_k\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right).\end{array}$$

Expressed in terms of the magnitude and direction of the mean wind vectors,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{w_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right.\right.\\ {\displaystyle }& {\displaystyle }\left.\left.-{\displaystyle \frac{\mathrm{2}\mathit{\rho }{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\right)\\ {\displaystyle }& {\displaystyle }\times \sum _{k=\mathrm{0}}^{\mathrm{\infty }}\left[{\mathit{ϵ}}_k\cos \left(k\mathit{\nu }\right)I_k\left({\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\sqrt{a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}}}\right)I_k\left({\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\sqrt{a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}}}\right)\right.\\ \text{(26)}& {\displaystyle }& {\displaystyle }\left.\times I_k\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right)\right],\end{array}$$

where

$$\begin{array}{ll}\text{(27)}& {\displaystyle }& {\displaystyle }\mathit{\nu }={\displaystyle }& {\displaystyle }{\tan }^{-\mathrm{1}}{\displaystyle \frac{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}\sin ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })}{(\mathrm{1}+{\mathit{\rho }}^{\mathrm{2}})\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })-\mathit{\rho }\left(\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}+\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\right)}}\end{array}$$

and

$$\begin{array}{ll}\text{(28)}& {\displaystyle }& {\displaystyle }\sqrt{a_i^{\mathrm{2}}+b_i^{\mathrm{2}}}={\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_i^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}}+{\mathit{\rho }}^{\mathrm{2}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{3}-i}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{3}-i}^{\mathrm{2}}}}-\mathit{\rho }{\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

When numerically evaluating the bivariate Rice distribution, it is convenient to transform the infinite series in Eqs. (25) and (26) into an integral which is then approximated numerically (Appendix A).

Note that p(w₁,w₂) depends on the relative orientation of the mean wind vectors and the correlation angle ψ only through the combination ${\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi }$. Because of this symmetry, the quantities ϕ₁−ϕ₂ and ψ cannot be determined individually from wind speed data alone. As a result, p(w₁,w₂) is determined by six parameters: $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}},\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi }$). For particular applications, it may be appropriate to fix either ϕ₁−ϕ₂ or ψ, allowing the other angle to be estimated from data. For example, when considering the temporal dependence structure of winds assumed to have stationary statistics, it can be assumed that ${\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}=\mathrm{0}$. The bivariate Rice distribution also has the discrete symmetry that it is invariant under the transformation ${\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi }\to -({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })$. Note that these symmetries are in addition to the invariance of the distribution of components (Eq. 5) to the rotation of the coordinate system ${\mathit{\theta }}_i\to {\mathit{\theta }}_i+\mathrm{\Delta }\mathit{\theta }$, $i=\mathrm{1},\mathrm{2}$, under which the angles ϕ₁−ϕ₂ and ψ are individually invariant.

Integrating over w₂ to obtain the marginal distribution for w₁ we obtain the univariate Rice distribution

$$\begin{array}{}\text{(29)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}p(w_{\mathrm{1}},w_{\mathrm{2}})\mathrm{d}{w}_{\mathrm{2}}={\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)I_{\mathrm{0}}\left({\displaystyle \frac{w_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right),\end{array}$$

with mean and variance

$$\begin{array}{}\text{(30)}& {\displaystyle }& {\displaystyle }\mathrm{mean}\left(w_{\mathrm{1}}\right)={\mathit{\sigma }}_{\mathrm{1}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}}{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\mathrm{1},-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)\text{(31)}& {\displaystyle }& {\displaystyle }\mathrm{var}\left(w_{\mathrm{1}}\right)=\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}{}_{\mathrm{1}}{F}_{\mathrm{1}}^{\mathrm{2}}\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\mathrm{1},-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}\right)\end{array}$$

(with equivalent expressions for w₂ obtained by integrating over w₁), where ${}_{\mathrm{1}}{F}_{\mathrm{1}}(\mathit{\alpha };\mathit{\beta };z)$ is the confluent hypergeometric function (Gradshteyn and Ryzhik, 2000). Equation (29) follows from Eq. (25) using the integral (Gradshteyn and Ryzhik, 2000)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}x{e}^{-a{x}^{\mathrm{2}}}{I}_k\left(bx\right)I_k\left(cx\right)\mathrm{d}x=\\ \text{(32)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}a}}\exp \left({\displaystyle \frac{b^{\mathrm{2}}+c^{\mathrm{2}}}{\mathrm{4}a}}\right)I_k\left({\displaystyle \frac{bc}{\mathrm{2}a}}\right),\end{array}$$

Neumann's theorem (Watson, 1922)

$$\begin{array}{}\text{(33)}& \sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\mathit{ϵ}}_k{I}_k\left(x\right)I_k\left(y\right)\cos k\mathit{\varphi }=I_{\mathrm{0}}\left(\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+\mathrm{2}xy\cos \mathit{\varphi }}\right),\end{array}$$

and the fact that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\cos \mathit{\nu }=\\ \text{(34)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{(\mathrm{1}+{\mathit{\rho }}^{\mathrm{2}})\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })-\mathit{\rho }\left(\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}+\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\right)}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}{)}^{\mathrm{2}}\sqrt{(a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}})\left(a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}}\right)}}}.\end{array}$$

Note that each wind speed marginal distribution depends only on the magnitude of the mean wind vector, while the joint distribution also depends on the angle between the two mean wind vectors. Furthermore, as ρ→0 only the first term contributes to the infinite series in Eq. (25), and the joint distribution reduces to the product of the marginals.

The joint moments of the bivariate Rice distribution can be evaluated using the Taylor series expansion:

$$\begin{array}{ll}\text{(35)}& {\displaystyle }& {\displaystyle }{I}_k\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right)={\displaystyle }& {\displaystyle }\sum _{j=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}}{j\mathrm{!}(j+k)\mathrm{!}}}{\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\right)}^{\mathrm{2}j+k}{\left({\displaystyle \frac{w_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{1}}}}\right)}^{\mathrm{2}j+k}{\left({\displaystyle \frac{w_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}}}\right)}^{\mathrm{2}j+k},\end{array}$$

which allows the double integral defining the moments to factorize as the products of individual integrals over w₁ and w₂ that can be evaluated using

$$\begin{array}{ll}\text{(36)}& {\displaystyle }& {\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{\mathit{\mu }}{e}^{-\mathit{\alpha }{x}^{\mathrm{2}}}{I}_k\left(\mathit{\beta }x\right)\mathrm{d}x={\displaystyle }& {\displaystyle }{\displaystyle \frac{{\mathit{\beta }}^k\mathrm{\Gamma }\left[(k+\mathit{\mu }+\mathrm{1})/\mathrm{2}\right]}{{\mathrm{2}}^{k+\mathrm{1}}{\mathit{\alpha }}^{(\mathit{\mu }+\mathit{\nu }+\mathrm{1})/\mathrm{2}}\mathrm{\Gamma }(k+\mathrm{1})}}{}_{\mathrm{1}}{F}_{\mathrm{1}}\left({\displaystyle \frac{k+\mathit{\mu }+\mathrm{1}}{\mathrm{2}}};k+\mathrm{1};{\displaystyle \frac{{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{4}\mathit{\alpha }}}\right).\end{array}$$

The resulting expression for the joint moments is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathsf{E}\left\{w_{\mathrm{1}}^m{w}_{\mathrm{2}}^n\right\}={\mathrm{2}}^{m/\mathrm{2}}{\mathrm{2}}^{n/\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}^m{\mathit{\sigma }}_{\mathrm{2}}^n(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}{)}^{\mathrm{1}+m/\mathrm{2}+n/\mathrm{2}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}\right.\right.\\ {\displaystyle }& {\displaystyle }\left.\left.-{\displaystyle \frac{\mathrm{2}\mathit{\rho }{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\right)\\ {\displaystyle }& {\displaystyle }\times \sum _{l=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^l\left[{\mathit{ϵ}}_k{\mathit{\rho }}^{\mathrm{2}l-k}\cos \left(k\mathit{\nu }\right){\displaystyle \frac{\mathrm{\Gamma }\left(l+\frac{m}{\mathrm{2}}+\mathrm{1}\right)\mathrm{\Gamma }\left(l+\frac{n}{\mathrm{2}}+\mathrm{1}\right)}{l\mathrm{!}(l-k)\mathrm{!}(k\mathrm{!}{)}^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\times {\left((\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{\sqrt{(a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}})(a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}})}}{\mathrm{2}}}\right)}^k\\ {\displaystyle }& {\displaystyle }\times {}_{\mathrm{1}}{F}_{\mathrm{1}}\left(l+{\displaystyle \frac{m}{\mathrm{2}}}+\mathrm{1};k+\mathrm{1};(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{(a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}})}{\mathrm{2}}}\right)\\ \text{(37)}& {\displaystyle }& {\displaystyle }\times \left.{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(l+{\displaystyle \frac{n}{\mathrm{2}}}+\mathrm{1};k+\mathrm{1};(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{(a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}})}{\mathrm{2}}}\right)\right]\end{array}$$

(Mendes and Yacoub, 2007). When the mean vector winds are equal to zero, only the k=0 terms contribute to this expression and Eq. (12) is recovered.

Figure 2 Comparison of the correlation coefficient corr(w₁,w₂) for bivariate Rice distributed variables (Eq. 38) with the approximate expression Eq. (40) for the parameter values $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}/{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{2}},{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{3},\mathrm{2},\mathrm{0})$ and $(\mathrm{0.5},\mathrm{4},\mathrm{0})$.

Defining the variables $V_i={\stackrel{\mathrm{‾}}{\mathcal{U}}}_i/\sqrt{\mathrm{2}}{\mathit{\sigma }}_i$, the correlation coefficient between w₁ and w₂ is given by

$$\begin{array}{ll}\text{(38)}& {\displaystyle }& {\displaystyle }\mathrm{corr}(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle }& {\displaystyle \frac{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}{)}^{\mathrm{2}}\exp \left(-\frac{V_{\mathrm{1}}^{\mathrm{2}}+V_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}{V}_{\mathrm{1}}{V}_{\mathrm{2}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}\right)G(V_{\mathrm{1}},V_{\mathrm{2}},\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })-\frac{\mathit{\pi }}{\mathrm{4}}{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{\mathrm{2}};\mathrm{1};-V_{\mathrm{1}}^{\mathrm{2}}\right){}_{\mathrm{1}}{F}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{\mathrm{2}};\mathrm{1};-V_{\mathrm{2}}^{\mathrm{2}}\right)}{{\left[\mathrm{1}+V_{\mathrm{1}}^{\mathrm{2}}-\frac{\mathit{\pi }}{\mathrm{4}}{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{\mathrm{2}};\mathrm{1};-V_{\mathrm{1}}^{\mathrm{2}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}{\left[\mathrm{1}+V_{\mathrm{2}}^{\mathrm{2}}-\frac{\mathit{\pi }}{\mathrm{4}}{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{\mathrm{2}};\mathrm{1};-V_{\mathrm{2}}^{\mathrm{2}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}}},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }G(V_{\mathrm{1}},V_{\mathrm{2}},\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=\\ {\displaystyle }& {\displaystyle }\sum _{l=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^l\left[{\mathit{ϵ}}_k{\mathit{\rho }}^{\mathrm{2}l-k}\cos \left(k\mathit{\nu }\right){\displaystyle \frac{{\mathrm{\Gamma }}^{\mathrm{2}}\left(l+\frac{\mathrm{3}}{\mathrm{2}}\right)}{l\mathrm{!}(l-k)\mathrm{!}(k\mathrm{!}{)}^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\times {\left((\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{\sqrt{(a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}})(a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}})}}{\mathrm{2}}}\right)}^k\\ {\displaystyle }& {\displaystyle }\times {}_{\mathrm{1}}{F}_{\mathrm{1}}\left(l+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}};k+\mathrm{1};(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{(a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}})}{\mathrm{2}}}\right)\\ \text{(39)}& {\displaystyle }& {\displaystyle }\times \left.{}_{\mathrm{1}}{F}_{\mathrm{1}}\left(l+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}};k+\mathrm{1};(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}){\displaystyle \frac{(a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}})}{\mathrm{2}}}\right)\right].\end{array}$$

The correlation coefficient corr(w₁,w₂) depends only on the four quantities $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}/{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{2}},{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi },\mathit{\rho })$. Monahan (2012b) considered the correlation structure of wind speeds using the approximation $\mathrm{corr}\left(w_{\mathrm{1}},w_{\mathrm{2}}\right)\simeq \mathrm{corr}\left(w_{\mathrm{1}}^{\mathrm{2}},w_{\mathrm{2}}^{\mathrm{2}}\right)$. The assumed covariance structure of the wind components then results in the approximate expression:

$$\begin{array}{ll}{\displaystyle }\mathrm{corr}\left(w_{\mathrm{1}},w_{\mathrm{2}}\right)& {\displaystyle }\simeq {\displaystyle \frac{{\mathit{\rho }}^{\mathrm{2}}+\left(\frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}+\frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}\right){\mathit{\mu }}_{\mathrm{1}}+\left(\frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}-\frac{{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}\right){\mathit{\mu }}_{\mathrm{2}}}{\sqrt{\left(\mathrm{1}+\frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}\right)\left(\mathrm{1}+\frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\right)}}}\\ \text{(40)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}^{\mathrm{2}}+\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}\cos \left({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi }\right)\mathit{\rho }}{\sqrt{\left(\mathrm{1}+\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}\right)\left(\mathrm{1}+\frac{{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\right)}}}.\end{array}$$

Plots of the correlation coefficient (Eq. 38) and the approximation (Eq. 40) as functions of ρ are shown in Fig. 2 for $(\stackrel{\mathrm{‾}}{{\mathcal{U}}_{\mathrm{1}}}/{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{2}},{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })$ = (3,2,0) and (1,5,0). Agreement between the exact and approximate values of the correlation coefficient is reasonably good in both cases, with the largest discrepancies generally occurring for larger absolute values of ρ. Note that negative wind speed correlation values are permitted by the bivariate Rice distribution.

Figure 3 As in Fig. 1 for the bivariate Rice distribution with $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}})=(\mathrm{6},\mathrm{4},\mathrm{2},\mathrm{5})$ and $(\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{0.8},\mathit{\pi }),(\mathrm{0},\mathrm{0}),$ and (0.8,0).

Examples of the joint Rice pdf (and the associated marginals) are presented in Fig. 3 for $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}})=(\mathrm{6},\mathrm{4},\mathrm{2},\mathrm{5})$ and $(\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{0.85},\mathit{\pi }),(\mathrm{0},\mathrm{0})$, and (0.85,0). By construction, the marginal distributions are the same in each panel. The distributions of both w₁ and w₂ are positively skewed, and take respective maxima at values of about σ₁ and just less than 2σ₂. For the different values of the dependence parameter ρ, the joint distributions have considerably different shapes. The joint distribution for $(\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{0.8},\mathit{\pi })$ describes (weakly) negatively correlated variables with a nonlinear dependence structure evident in ridges of enhanced probability extending to the left and right upward from the probability maximum. For $(\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{0},\mathrm{0})$, probability contours are concentrated towards smaller values of w₁ and w₂ (as is the case for the marginal distributions). Finally, w₁ and w₂ are evidently positively correlated for $(\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{0.8},\mathrm{0})$, with a slight curvature in the shape of the distribution indicating the existence of some nonlinear dependence.

Although the bivariate Rice distribution differs from the bivariate Rayleigh distribution only by allowing for nonzero mean wind vector components, the resulting expressions for the joint pdf (Eq. 25) and the moments (Eq. 37) are much more complicated for the bivariate Rice than the bivariate Rayleigh. Furthermore, while the univariate Rice distribution is a convenient model for the pdf of wind speed, observed winds show clear deviation from Ricean behaviour (e.g. Monahan, 2006, 2007). We will therefore consider another model of the bivariate wind speed distribution with Weibull marginals, which turns out to result in simpler mathematical expressions (at the cost of a more artificial derivation than that of the bivariate Rice distribution).

2.3 Bivariate Weibull distribution

As in Sagias and Karagiannidis (2005) and Yacoub et al. (2005), we obtain the bivariate Weibull distribution from the bivariate Rayleigh distribution through separate power law transformations of w₁ and w₂. The pdf of a Weibull distributed variable is

$$\begin{array}{}\text{(41)}& p\left(x\right)={\displaystyle \frac{b}{a}}{\left({\displaystyle \frac{x}{a}}\right)}^{b-\mathrm{1}}\exp \left(-{\left[{\displaystyle \frac{x}{a}}\right]}^b\right),\end{array}$$

with moments

$$\begin{array}{}\text{(42)}& \mathsf{E}\left\{x^m\right\}=a^m\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{m}{b}}\right),\end{array}$$

where a and b are denoted the scale and shape parameters respectively. The Rayleigh distribution is a special case of the Weibull distribution with $a=\sqrt{\mathrm{2}}\mathit{\sigma }$ and b=2. Weibull distributed variables remain Weibull under a power law transformation, with suitably modified scale and shape parameters: if x is Weibull with scale parameter a and shape parameter b, x^k will be Weibull with scale parameter a^k and shape parameter b∕k. A joint wind speed distribution with Weibull marginal distributions can therefore be constructed from a joint Rayleigh distribution using the appropriate power law and scale transformations.

If we start with (x₁,x₂) as bivariate Rayleigh distributed with ${\mathit{\sigma }}_i=\mathrm{1}/\sqrt{\mathrm{2}}$, $i=\mathrm{1},\mathrm{2}$, we obtain marginal Weibull distributions with specified scale and shape parameters through the transformations

$$\begin{array}{}\text{(43)}& w_i=a_i{x}_i^{\mathrm{2}/b_i}.\end{array}$$

The joint pdfs transform as

$$\begin{array}{}\text{(44)}& p(w_{\mathrm{1}},w_{\mathrm{2}})={\left|\begin{array}{cc}{\partial }_{x_{\mathrm{1}}}{w}_{\mathrm{1}}& {\partial }_{x_{\mathrm{2}}}{w}_{\mathrm{1}}\\ {\partial }_{x_{\mathrm{1}}}{w}_{\mathrm{2}}& {\partial }_{x_{\mathrm{2}}}{w}_{\mathrm{2}}\end{array}\right|}^{-\mathrm{1}}p(x_{\mathrm{1}},x_{\mathrm{2}})\end{array}$$

and so we obtain the bivariate Weibull distribution

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{b_{\mathrm{1}}{b}_{\mathrm{2}}}{a_{\mathrm{1}}{a}_{\mathrm{2}}}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}-\mathrm{1}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}}+{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}}\right]\right)\\ \text{(45)}& {\displaystyle }& {\displaystyle }\times I_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right).\end{array}$$

An analogous approach to constructing bivariate Weibull distributions through nonlinear transformations of a bivariate Gaussian was followed in Villanueva et al. (2013); the resulting expressions are considerably more complicated than those considered here.

Evidently, p(w₁,w₂) factorizes into the product of the marginal distributions as ρ→0. As ρ→1, we can use Eq. (10) to make the approximation

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(w_{\mathrm{1}},w_{\mathrm{2}})\simeq {\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{b_{\mathrm{1}}{b}_{\mathrm{2}}}{a_{\mathrm{1}}{a}_{\mathrm{2}}}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}-\mathrm{1}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\left[{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}}+{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}}\right]\right)\\ {\displaystyle }& {\displaystyle }\times {\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{2}}\\ {\displaystyle }& {\displaystyle }\times \exp \left({\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\simeq \left[{\displaystyle \frac{b_{\mathrm{1}}}{a_{\mathrm{1}}}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}-\mathrm{1}}\exp \left(-{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}}\right)\right]\\ {\displaystyle }& {\displaystyle }\times \left[{\displaystyle \frac{b_{\mathrm{2}}}{\mathrm{2}{a}_{\mathrm{2}}}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{-b_{\mathrm{1}}/\mathrm{4}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{\mathrm{3}{b}_{\mathrm{2}}/\mathrm{4}-\mathrm{1}}\right]\\ {\displaystyle }& {\displaystyle }\times \mathit{\delta }\left({\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}-{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right)\\ \text{(46)}& {\displaystyle }& {\displaystyle }=p\left(w_{\mathrm{1}}\right)\mathit{\delta }\left(w_{\mathrm{2}}-a_{\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/b_{\mathrm{2}}}\right),\end{array}$$

where the last equality follows from the fact that

$$\begin{array}{}\text{(47)}& \mathit{\delta }\left(f\left(x\right)-\mathit{\alpha }\right)={\displaystyle \frac{\mathrm{1}}{f^{\prime }\left(\mathit{\alpha }\right)}}\mathit{\delta }\left(x-f^{-\mathrm{1}}\left(\mathit{\alpha }\right)\right).\end{array}$$

As expected, w₁ and w₂ are completely dependent in the limit that ρ→1 (although they are not perfectly correlated if b₁≠b₂ as the functional relationship

$$\begin{array}{}\text{(48)}& w_{\mathrm{2}}=a_{\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/b_{\mathrm{2}}}\end{array}$$

between the two variables will be nonlinear).

The relatively simple form of the bivariate Weibull distribution permits a relatively simple expression for the conditional distribution

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p\left(w_{\mathrm{2}}\right|w_{\mathrm{1}})={\displaystyle \frac{p(w_{\mathrm{1}},w_{\mathrm{2}})}{p\left(w_{\mathrm{1}}\right)}}\\ {\displaystyle }& {\displaystyle }=\left\{{\displaystyle \frac{b_{\mathrm{2}}}{a_{\mathrm{2}}}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}-\mathrm{1}}\exp \left(-{\left[{\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right]}^{b_{\mathrm{2}}}\right)\right\}\\ {\displaystyle }& {\displaystyle }\times \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}}+{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}}\right]\right)\right.\\ \text{(49)}& {\displaystyle }& {\displaystyle }\left.\times I_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{2}\mathit{\rho }}{(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right)\right\}.\end{array}$$

The factor in the second set of braces characterizes how conditioning on the value of w₁ changes the distribution of w₂ from its marginal distribution (and corresponds to a copula density; e.g. Schlözel and Friederichs, 2008). Note that for w₁ sufficiently large we can write

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p\left(w_{\mathrm{2}}\right|w_{\mathrm{1}})\simeq {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{4}\mathit{\pi }\mathit{\rho }(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}}{\displaystyle \frac{b_{\mathrm{2}}}{a_{\mathrm{2}}}}{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{\mathrm{3}{b}_{\mathrm{2}}/\mathrm{4}-\mathrm{1}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{-b_{\mathrm{1}}/\mathrm{4}}\\ \text{(50)}& {\displaystyle }& {\displaystyle }\times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\left[\mathit{\rho }{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}/\mathrm{2}}-{\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}/\mathrm{2}}\right]}^{\mathrm{2}}\right).\end{array}$$

For ρ not too close to zero, the conditional distribution for large w₁ is concentrated around the nonlinear regression curve

$$\begin{array}{}\text{(51)}& {\left({\displaystyle \frac{w_{\mathrm{2}}}{a_{\mathrm{2}}}}\right)}^{b_{\mathrm{2}}}={\mathit{\rho }}^{\mathrm{2}}{\left({\displaystyle \frac{w_{\mathrm{1}}}{a_{\mathrm{1}}}}\right)}^{b_{\mathrm{1}}}.\end{array}$$

Figure 4 As in Fig. 1 for the bivariate Weibull distribution with $(a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}})=(\mathrm{4},\mathrm{1.5},\mathrm{5},\mathrm{7})$ and $\mathit{\rho }=\mathrm{0},\mathrm{0.7}$, and 0.95 and the speeds w₁,w₂ scaled respectively by the Weibull scale parameters a₁, a₂.

Computing the moments, we obtain

$$\begin{array}{ll}\text{(52)}& {\displaystyle }& {\displaystyle }\mathsf{E}\left\{w_{\mathrm{1}}^m{w}_{\mathrm{2}}^n\right\}={\displaystyle }& {\displaystyle }{a}_{\mathrm{1}}^m{a}_{\mathrm{2}}^n\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{m}{b_{\mathrm{1}}}}\right)\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{n}{b_{\mathrm{2}}}}\right){}_{\mathrm{2}}{F}_{\mathrm{1}}\left(-{\displaystyle \frac{m}{b_{\mathrm{1}}}},-{\displaystyle \frac{n}{b_{\mathrm{2}}}},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}}\right).\end{array}$$

(Sagias and Karagiannidis, 2005; Yacoub et al., 2005). The correlation coefficient is then given by

$$\begin{array}{ll}\text{(53)}& {\displaystyle }& {\displaystyle }\mathrm{corr}(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{1}+\frac{\mathrm{1}}{b_{\mathrm{1}}}\right)\mathrm{\Gamma }\left(\mathrm{1}+\frac{\mathrm{1}}{b_{\mathrm{2}}}\right)\left[{}_{\mathrm{2}}{F}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{b_{\mathrm{1}}},-\frac{\mathrm{1}}{b_{\mathrm{2}}},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}}\right)-\mathrm{1}\right]}{\sqrt{\left[\mathrm{\Gamma }\left(\mathrm{1}+\frac{\mathrm{2}}{b_{\mathrm{1}}}\right)-{\mathrm{\Gamma }}^{\mathrm{2}}\left(\mathrm{1}+\frac{\mathrm{1}}{b_{\mathrm{1}}}\right)\right]\left[\mathrm{\Gamma }\left(\mathrm{1}+\frac{\mathrm{2}}{b_{\mathrm{2}}}\right)-{\mathrm{\Gamma }}^{\mathrm{2}}\left(\mathrm{1}+\frac{\mathrm{1}}{b_{\mathrm{2}}}\right)\right]}}}.\end{array}$$

For $b_{\mathrm{1}},b_{\mathrm{2}}>\mathrm{1}$ (the relevant range of shape parameters for wind speeds), ${}_{\mathrm{2}}{F}_{\mathrm{1}}(-\mathrm{1}/b_{\mathrm{1}},-\mathrm{1}/b_{\mathrm{2}},\mathrm{1};{\mathit{\rho }}^{\mathrm{2}})$ is an increasing function of ρ² with ${}_{\mathrm{2}}{F}_{\mathrm{1}}(-\mathrm{1}/b_{\mathrm{1}},-\mathrm{1}/b_{\mathrm{2}},\mathrm{1};\mathrm{0})=\mathrm{1}$. Therefore, the bivariate Weibull distribution is unable to represent situations in which the wind speeds are negatively correlated.

Examples of the bivariate Weibull distribution for $(a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}})=(\mathrm{4},\mathrm{1.5},\mathrm{5},\mathrm{7})$ and values of $\mathit{\rho }=\mathrm{0},\mathrm{0.7}$, and 0.95 are shown in Fig. 4. Again, the marginal distributions in all three cases are the same by construction. The distribution of w₁ is positively skewed with a maximum near a value of w₁=0.5a₁, while that of w₂ is negatively skewed with a maximum near w₂=a₂. For ρ=0 the joint distribution is simply the product of the marginals. As ρ increases, w₁ and w₂ become positively correlated - although the correlation is weak even for ρ=0.7 for this set of parameter values. At the value of ρ=0.95, while the correlation of the two variables is only moderate, a strong nonlinear dependence is evident in the concentration of the distribution around the curve given by Eq. (48).

Table 1 Maximum likelihood parameter estimates for the wind speed data shown in Figs. 5, 6, 7, and 8.

3 Fits of bivariate Rice and Weibull distributions to observed wind speeds

Many wind datasets from different locations are available, and it is impracticable to consider joint distributions of wind speeds from even a small fraction of these. In this section, we will consider examples of the joint distribution of wind speeds using data from a representative range of settings. Bivariate distributions of wind speeds at both different locations in space and different points in time will be considered. The sampling of the wind speeds considered will be temporal (that is, individual samples will correspond to a specific time for spatial joint pdfs and a specific pair of times for temporal joint pdfs). Best-fit values of the parameters of the bivariate Weibull and Rice distributions we present were obtained numerically as maximum likelihood estimates (Table 1), with ϕ₁−ϕ₂ set to zero. Goodness-of-fit of the distributions was assessed using a statistical test described in Appendix B. In order to distinguish how well the parametric joint distributions model the marginal distributions from how well they represent dependence between variables, the goodness-of-fit analyses were repeated for each pair of time series with the values of one of the pair shuffled in time. This shuffling destroys the dependence structure without affecting the distributions of the marginals. Use of a bivariate analysis rather than separate univariate goodness-of-fit tests for the marginals allows direct comparisons of p-values, as exactly the same test is used for the original and shuffled data.

Figure 5 Joint distributions of 500 hPa DJF 00Z wind speeds at four different pairs of latitudes along 216^∘ W. Wind speeds at the two latitudes quoted in the column headings are respectively denoted w₁ and w₂. (a–d) Scatterplots of wind speed data. (e–h) Maximum likelihood bivariate Weibull pdfs (white contours) and kernel density estimates of the observed joint pdf (colours). The p-value of a goodness-of-fit test with the null hypothesis that the observed wind speed data are drawn from the corresponding best-fit bivariate Weibull distribution is given. The values in brackets are the p-values obtained when the dependence structure is eliminated by shuffling the values of w₂ in time. (i–l) As in the middle row, for the best-fit bivariate Rice distribution. Values of the best-fit model parameters are given in Table 1.

3.1 Wind speeds at 500 hPa

We first consider the joint distribution of 00Z December, January, and February 500 hPa wind speeds from 1979 to 2014. These data were taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim Reanalysis (Dee et al., 2011), subsampled to every second day to minimize the effect of serial dependence on the goodness-of-fit test (Monahan, 2012b). The wind speed data were computed as the magnitude of zonal and meridional components.

The joint distributions of wind speeds at four pairs of latitudes along 216^∘ W are presented in Fig. 5. A moderately strong negative correlation ($r=-\mathrm{0.55}$) is evident between wind speeds at (39^∘ S, 216^∘ W) and (54.75^∘ S, 216^∘ W) (Fig. 5a). Because it is unable to model a negative correlation between wind speeds, the best-fit bivariate Weibull pdf differs substantially from the distribution of the observed winds (Fig. 5e). The goodness-of-fit test correspondingly rejects the null hypothesis that the observations are drawn from this distribution (p=0). In contrast, the bivariate Rice distribution provides a reasonable model of the joint distribution of wind speeds at these two locations (Fig. 5i) and the goodness-of-fit test provides no evidence that these data are statistically incompatible with this distribution (p=0.31). For wind speeds at these two locations, the bivariate Rice distribution is evidently a better model than the bivariate Weibull distribution.

In contrast, for wind speeds at (12^∘ N, 216^∘ W) and (15.75^∘ N, 216^∘ W), the null hypotheses of being drawn from either the bivariate Rice or Weibull distributions are rejected at the 95 % significance level for both distributions. These wind speeds are weakly correlated (r=0.37) but show evidence of nonlinear dependence. The joint pdf of the observed speeds is characterized by two ridges of high probability extending to the right, upward and downward away from the region of maximum probability (Fig. 5b, f, j). These ridges are not captured by either of the best-fit bivariate Weibull or Rice distributions, although a hint of this structure is evident in the Rice distribution. While the fit of one or both of these parametric distributions to these observed wind speed data may be sufficiently good for practical applications, nevertheless we can confidently exclude the possibility that these data are drawn from either distribution.

The wind speeds at (3^∘ S, 216^∘ W) and (0.75^∘ N, 216^∘ W) are correlated (r=0.68) and their scatter clusters around a straight line extending away from the origin (Fig. 5c). Both the best-fit bivariate Weibull and Rice distributions appear to the eye to be good fits to the data (Fig. 5g, k), and in neither case can the null hypothesis be rejected that the data are drawn from these distributions. Only small differences exist between the two best-fit distributions for these data.

Finally, the wind speeds at (15^∘ S, 216^∘ W) and (45^∘ S, 216^∘ W) are uncorrelated (r=0.06) and fit sufficiently well by both the bivariate Weibull and Rice distributions (Fig. 5d, h, l) that in neither case is the null hypothesis rejected. As in the previous example, the bivariate best-fit Rice and Weibull distributions are essentially indistinguishable for these data.

Figure 6 The temporal dependence structure of 500 hPa DJF 00Z wind speeds at (39^∘ S, 216^∘ W) for lags s = 0.25 to 4 days, with w₁=w(t_n) and $w_{\mathrm{2}}=w\left(t_{n+s}\right)$. (a–c) Scatterplots of wind speeds separated by 0.5 days (first column), 1.5 days (second column), and 3 days (third column). (e–g) Maximum likelihood bivariate Weibull distribution (white contours) and kernel density estimate of the joint pdf of the lagged data. The p-value of a goodness-of-fit test of the bivariate Weibull fit is quoted in white. The values in brackets are the p-values obtained when the dependence structure is eliminated by shuffling the values of w₂ in time. (i–k) As in panels (e–g) for the bivariate Rice distribution. Panels (d, h, l) show results at a range of different time lags. (d) p-values of bivariate goodness-of-fit tests for bivariate Weibull (solid black) and bivariate Rice (solid red). The dashed lines correspond to the p-values for shuffled w₂. The thin black curve is the 0.95 significance level. (h) Best-fit estimate of the parameter ρ of the bivariate Weibull distribution. (l) Best-fit estimates of ρ (black) and cos ϕ (blue) for the Rice distribution. Values of the best-fit model parameters are given in Table 1.

Considering the spatial correlation structure of these 500 hPa winds, we find cases in which one distribution (the Rice) is evidently a better fit to the data than the other (the Weibull), in which neither distribution provides a statistically significant fit to the data, and in which both distributions fit the data equally well. p-values from the analyses with temporally shuffled data, quoted in parentheses in Fig. 5, show that for none of the four pairs of wind speeds considered can the null hypotheses of univariate Weibull or univariate Rice distributions for the pair of marginals be rejected. This result indicates that for these data the rejection of the full bivariate distributions occurs because of a failure to adequately represent the dependence structure between the two random variables.

The temporal dependence structure of the wind speed at (39^∘ S, 216^∘ W) is illustrated in Fig. 6. As with the previous calculations, the pairs of lagged wind speeds $\left(w\right(t_n),w(t_{n+s}\left)\right)$ were subsampled to 2-day resolution to minimize the effect of serial dependence on the results of the goodness-of-fit tests. As the lag increases, the value of the dependence parameter ρ decreases as expected for both the Weibull and Rice distributions (Fig. 6h, l). For most lags, the null hypothesis of the Weibull distribution as a model for the joint distribution is rejected (p<0.05; Fig. 6d). The rejection of the null hypothesis of a bivariate Weibull distribution is most robust for lags shorter than 3 days. In contrast, the null hypothesis of a bivariate Rice distribution is rejected less often than it is not – although p<0.05 for more than one-third of the lags (Fig. 6d). Inspection of the example distributions shown demonstrates that the bivariate Weibull distributions are broader around their principal axis for small to intermediate wind speeds in a way that is not consistent with observations (Fig. 6e–g). Such structures are not seen in the best-fit Rice distributions (Fig. 6i–k). Note that for lags of 0.5 and 1.5 days the observed distributions suggest a flaring out of the joint distribution for large wind speeds that is accounted for by neither the bivariate Weibull nor Rice distributions. There is good evidence that these data were not drawn from a bivariate Weibull distribution, and the evidence that they are drawn from a bivariate Rice distribution is not strong. p-values of fits obtained after shuffling of the lagged wind speed time series (in brackets in Fig. 6e–g, i–k and dashed lines in Fig. 6d) are generally larger than those of the unshuffled data: with a few exceptions (e.g. a lag of 1.75 days), the rejection of the null hypothesis of the full dataset being drawn from either of the parametric distributions considered is not associated with a failure to fit the marginals.

Figure 7 Joint distribution of JAS wind speeds at 10 and 200 m measured at Cabauw, NL. (a–d) Wind speed scatterplots. (e–h) Kernel density estimate of the joint pdf of wind speeds (colour) and maximum likelihood bivariate Weibull distribution (contours). The p-value of a goodness-of-fit test of the bivariate Weibull distribution is quoted in white. The values in brackets are the p-values obtained when the dependence structure is eliminated by shuffling the values of w₂₀₀ in time. (i–l) As in the middle row, but for the best-fit bivariate Rice distribution (contours). (a, e, i) All nighttime data (20:00–05:00 UTC). (b, f, j) Nighttime data conditioned on being in regime R₂ (very stable boundary layer). (c, g, k) Nighttime data conditioned on being in regime R₁ (weakly stable boundary layer). (d, h, l) All daytime data (08:00–16:00 UTC). Values of the best-fit model parameters are given in Table 1.

3.2 Wind speeds over land at 10 m and 200 m

We next consider wind speeds at altitudes of 10 and 200 m measured from a 213 m tower in Cabauw, Netherlands (51.971^∘ N, 4.927^∘ E), maintained by the Cabauw Experimental Site for Atmospheric Research (CESAR; van Ulden and Wieringa, 1996) with 10 min resolution from 1 January 2001 through 31 December 2012. We will focus on data from July, August, and September (JAS), separated into daytime (08:00–16:00 UTC) and nighttime (20:00–05:00 UTC) periods. These data were subsampled in time to account for serial dependence. Only every 50th point was used in the following analysis. A small number of zero wind speed values were removed from the dataset.

Monahan et al. (2011, 2015) demonstrated the existence of two distinct regimes of the nocturnal boundary layer in these data, corresponding to the very and weakly stratified boundary layers (vSBL and wSBL; e.g. Mahrt, 2014). These regimes, denoted respectively R₁ and R₂, were separated in Monahan et al. (2015) using a two-state HMM. Conditioning the data on the HMM state, the scatterplot of wind speeds at 10 and 200 m separates into two distinct populations (Fig. 7a–c). The scatter of wind speeds at these two altitudes shows no evident regime structure during the day (Fig. 7d). In all cases, the wind speeds at the two altitudes are highly correlated.

Maximum likelihood estimates of the bivariate Weibull and Rice distributions for the full nighttime data show evident disagreement between the scatter of the data and the best-fit distributions (Fig. 7e, i). Goodness-of-fit tests for both distributions were rejected (with p=0 in both cases). When conditioned on being in either regime R₁ or R₂, both the bivariate Rice and Weibull distributions result in much better representations of the nighttime data (Fig. 7f, g, j, k). In all cases the p-values exceed 0.05, so the fits cannot be rejected at the 95 % significance level. While neither the bivariate Rice nor Weibull distributions is a good probability model for the full joint distribution of wind speeds at these two altitudes, both are reasonable models for the distributions conditioned by regime occupation. Finally, the fits of neither the bivariate Rice nor Weibull distributions to the daytime data are statistically significant at the 95 % significance level. In particular, the joint Rice distribution is too broad (relative to observations) for small values of w₁₀ and w₂₀₀ and neither distribution is broad enough at higher wind speed values. p-values of fits with w₂₀₀ shuffled in time are all larger than those of fits to the original data; only for the Rice fit to the full nighttime data is the shuffled p-value below 0.05. As seen in the previous data considered, there is no systematic evidence that the failure of the joint Weibull or Rice distributions to model the joint distributions of the Cabauw data results from a failure to model the marginals.

Figure 8 Joint distribution of DJF QuikSCAT wind speeds at (6.5^∘ S, 135^∘ W), denoted w₁, and westward along a zonal transect to 165^∘ W (denoted w₂). (a–c) Scatterplots of wind speed at (6.5^∘ S, 135^∘ W) and at the points (6.5^∘ S, 162^∘ W), (6.5^∘ S, 152^∘ W), and (6.5^∘ S, 142^∘ W). (e–g) Kernel density estimate of the joint pdf (colour) as well as the maximum likelihood bivariate Weibull distribution (white contours). The p-value of a goodness-of-fit test of the bivariate Weibull distribution is quoted in white. The values in brackets are the p-values obtained when the dependence structure is eliminated by shuffling the values of w₂ in time. (i–k) As in panels (e–g)  but for the bivariate Rice distribution. (d) p-values of the bivariate goodness-of-fit tests for the wind speeds along the transect, for the bivariate Weibull distribution (solid black) and for the bivariate Rice distribution (solid red). The dashed lines correspond to the p-values for shuffled w₂. (h) Estimate of the parameter ρ from the best-fit bivariate Weibull distribution along the transect. (l) Estimate of the parameters ρ (black) and cos ψ (blue) from the best-fit bivariate Rice distribution along the transect. Values of the best-fit model parameters are given in Table 1.

3.3 Sea surface wind speeds

Twice-daily December, January, and February level 3.0 gridded SeaWinds scatterometer equivalent neutral 10 m wind speeds between 60^∘ S and 60^∘ N at a resolution of 0.25^∘ × 0.25^∘ from the National Aeronautics and Space Administration (NASA) Quick Scatterometer (QuickSCAT; Perry, 2001) are available from December 1999 through February 2008. Data flagged as having possibly been corrupted by rain were excluded from the following analysis. Although the data are nominally twice-daily, it is often the case that data for either the ascending or descending pass of the satellite are missing. The maximum likelihood parameter estimates of bivariate wind speed distributions and goodness-of-fit tests were carried out using every third non-missing data point in order to minimize the effect of serial dependence. For the goodness-of-fit tests M=10 quantiles were used because of the relatively small sample sizes. Because of the near-polar orbit of the satellite results, observations of wind speed at different locations are not simultaneous. The joint distributions we consider therefore combine dependence in both space and time.

Joint distributions of wind speed at (6.5^∘ S, 135^∘ W) with speeds at points along a zonal transect to 165^∘ W (in increments of 1^∘) were estimated. As the distance between the two positions increases, there is a decreasing trend in the best-fit bivariate Weibull dependence parameter ρ (Fig. 8h) with some small fluctuations likely due to sampling variability. The same is not true for the best-fit bivariate Rice dependence parameter ρ, which fluctuates wildly (Fig. 8l). Large fluctuations also seen in the best-fit value of cos ψ are clearly correlated with those of ρ: where one parameter is anomalously large (relative to the spatial trend), the other is anomalously small. Of the 30 pairs of points considered, the null hypothesis of a bivariate Rice distribution is rejected (at the 95 % significance level) at only one (Fig. 8d). In contrast, p<0.05 for the bivariate Weibull distribution at several longitudes, particularly close to the base point at (6.5^∘ S, 135^∘ W). Inspection of the best-fit bivariate Weibull pdfs (Fig. 8e–g) shows that they are broader for smaller wind speeds than for larger values, a feature not evident in the best-fit bivariate Rice distributions (Fig. 8i–k) or the scatter of data (Fig. 8a–c). From these results, we find only equivocal evidence that the pairs of wind speed data along this zonal transect are drawn from a bivariate Weibull distribution and no strong evidence to reject the null hypothesis that they are drawn from a bivariate Rice distribution. Again, the p-values of fits to temporally shuffled wind speeds are generally similar to or larger than those of fits to the full distribution. The bivariate Rice distribution with the wind speed at 142^∘ W (Fig. 8k) illustrates a rare example in which the parametric fit to the original data passed the goodness-of-fit test at the 5 % significance level while the fit to the shuffled data did not; it is evident from Fig. 8d that this situation is not common. The surface wind vector components in the tropics are known to be non-Gaussian (e.g. Monahan, 2007), so we have a priori reasons to believe the joint distribution should not be Ricean. The fact that the data do not generally allow for a rejection of the null hypothesis that the winds are bivariate Rice (for either the original or shuffled datasets) is likely a consequence of the relatively small sample size.

Figure 9 (a–c) Estimates of ρ and cos ψ from 50 realizations of bivariate Ricean variables with $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}},\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{7.3},\mathrm{1.5},\mathrm{6.9},\mathrm{1.5},\mathrm{0.4},\mathrm{0})$ for each of the sample sizes $N=\mathrm{250},\mathrm{1250}$, and 9000. The open circles correspond to estimates with $\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })\ge \mathrm{0.99}$, while the stars are estimates with $\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })<\mathrm{0.99}$. (d–f) Contours of bivariate Rice pdfs corresponding to 10 of the 50 best-fit parameter estimates (randomly chosen). The contour values are the same for all pdfs within each subplot. (g–i) Scatterplot of the sample correlation coefficient between the two Ricean variables and the correlation coefficient given by the approximate expression Eq. (40). The 1 : 1 line is given in solid black. The open circles and stars are as in the upper row.

The large variations in best-fit estimates of ρ and cos ψ for the bivariate Rice distribution result from the fact that for some parameter values the distribution is only weakly sensitive to simultaneous changes in these parameters: increases in ρ can be offset by decreases in cos ψ with only small changes to the joint distribution. To demonstrate this weak sensitivity, 50 realizations of bivariate Ricean variables with $({\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{2}},\mathit{\rho },{\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })=(\mathrm{7.3},\mathrm{1.5},\mathrm{6.9},\mathrm{1.5},\mathrm{0.5},\mathrm{0})$ were generated for each of the sample sizes of N=250, 1500, and 9000. Maximum likelihood estimates of these parameters obtained from these realizations demonstrate that for the smaller samples ρ and cos ψ show strong and correlated sampling variability, with large increases in ρ combined with large decreases in cos ψ (Fig. 9a). As expected, these sampling fluctuations become smaller as the sample size increases(Fig. 9b, c). Despite the large variation of ρ and cos ψ for small to intermediate sample sizes, there is relatively little variation in the structure of the corresponding bivariate Rice distributions (Fig. 9d–f).

An indication of why increases in ρ should counterbalance decreases in cos ψ with only small effects on the joint distribution is given by the approximate expression for the correlation coefficient, Eq. (40). The value of this approximation is unaffected by changes in ρ and ψ that leave the numerator invariant. The compensation between sampling variations in ρ and cos ψ is evident the fact that corr(w₁,w₂) given by Eq. (40) is an excellent approximation to the sample correlation coefficient even for estimates of ρ and cos ϕ which are far away from the population values (Fig. 9g–i). Note that there is no evident relationship between sampling fluctuations in corr(w₁,w₂) and those of ρ and ψ: the range of sample correlation values for (ρ,cos ψ) near the population values of (0.5,1) (open circles) is the same as that for values of (ρ,cos ψ) far from these values (stars). In these parameter ranges, the dependence between w₁ and w₂ constrains ρ and ψ not individually, but together – over large ranges of values for sufficiently small sample sizes.

4 Conclusions

This study has considered two idealized probability models for the joint distribution of wind speeds, both derived from models for the joint distribution of the horizontal wind components. The first, the bivariate Rice distribution, follows from assuming that the wind vector components are bivariate Gaussian with an idealized covariance structure. The second, the bivariate Weibull distribution, arises from nonlinear transformations of variables with a bivariate Rice distribution in the limit that the mean vector winds vanish (the bivariate Rayleigh distribution). While the bivariate Rice distribution has the advantage of being more flexible and naturally related to a simplified model for the joint distribution of the wind components, the bivariate Weibull distribution is mathematically much simpler and easier to work with. Through consideration of a range of joint distributions of observed wind speeds (over land and over the ocean; at the surface and aloft; in space and in time) the bivariate Rice distribution was shown to generally model the observations better than the bivariate Weibull distribution. However, in many circumstances the differences between the two distributions are small and the convenience of the bivariate Weibull distribution relative to the bivariate Rice distribution is a factor which may motivate its use.

The fact that the bivariate Rice distribution is easier to work with, but less flexible, than the bivariate Weibull distribution is evident from inspection of their analytic forms and the relative number of parameters to fit (five vs. six). If the bivariate Weibull distribution was generically appropriate for modelling the bivariate wind speed distribution, there would be no need to consider more complicated models such as the bivariate Rice. This study provides an empirical assessment of the relative practical utility of the two models, trading off the ability to model more general dependence structures (e.g. negatively correlated speeds) against model simplicity. Neither the univariate nor the bivariate Weibull or Rice distributions are expected to represent the true distributions of wind speeds (e.g. Carta et al., 2009). The results of this analysis characterize the practical utility of these models, rather than making a claim to their “truth”. It is noteworthy that for the data considered in this study, the failure of either the bivariate Rice or Weibull distributions to adequately fit the joint distribution of wind speeds (at a significance level of 5 %) is not generally associated with a corresponding failure of the parametric distribution to model the marginals. An interesting direction of future study would be consideration of other parametric models for the joint distribution of non-negative quantities, such as the α-μ distribution (Yacoub, 2007), copula-based models (e.g. Schlözel and Friederichs, 2008), or distributions obtained through nonlinear transformations of multivariate Gaussians (e.g. Brown et al., 1984).

Many of the assumptions that have been made regarding the distribution of the wind components are known not to hold in various settings. For example, the vector wind components are generally not Gaussian, either aloft or at the surface (e.g. Monahan, 2007; Luxford and Woollings, 2012; Perron and Sura, 2013), and fluctuations will not generally be isotropic (especially over land; cf. Mao and Monahan, 2017). Furthermore, when used to model temporal dependence the assumed correlation structure cannot account for the anisotropy in autocorrelation of orthogonal components in either space (e.g. Buell, 1960) or time (e.g. Monahan, 2012b). Relaxing the assumptions regarding isotropy of correlation structure results in expressions for the joint speed distributions involving integrals over angle which are not analytically tractable.

While it is possible to relax the assumption of Gaussian components for univariate speed distribution (e.g. Monahan, 2007; Drobinski et al., 2015), extending this analysis to the bivariate case would involve specifying a non-Gaussian dependence structure for the components. At present, there is no physically based model for such dependence. Without any such physical justification, the only option is an empirical investigation of the ability of different copula models to represent observed joint wind speed distributions. It is unlikely that a copula-based model for dependence of components will admit analytically tractable expressions for joint speed distributions. A copula-based analysis of either the components or the speeds directly is also likely necessary for modelling extreme wind speeds (either large percentiles, peaks over threshold, or block maxima), as the tails of the bivariate Rice and Weibull distributions may not be adequate for this task. Finally, extending the approach used in this study to obtain explicit closed-form results for the bivariate wind speed distribution to a higher-dimensional multivariate setting – of wind speeds alone, or of a mixture of wind speeds and other meteorological quantities – will be analytically intractable for any except the simplest (and likely unrealistic) covariance structures. It may not be practical to extend the program of obtaining closed-form expressions for joint speed distributions from models of the component distributions much further beyond the bivariate Rice and Weibull speed distributions considered in this study.

Ultimately, it would be best for models of the joint distribution of wind speeds to arise from physically based (if still idealized) models, as has been done for the univariate case in Monahan (2006) and Monahan et al. (2011). The development of such models represents an interesting direction of future study.

Appendix A: Numerical computation of bivariate Rice pdf

Equation (25) is difficult to evaluate numerically when the arguments of the Bessel functions become large. We have found that a computationally more stable result is obtained when this equation is expressed in the form

$$\begin{array}{}\text{(A1)}& p(w_{\mathrm{1}},w_{\mathrm{2}})={\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\exp f(w_{\mathrm{1}},w_{\mathrm{2}},\mathit{\lambda })\mathrm{d}\mathit{\lambda },\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }f(w_{\mathrm{1}},w_{\mathrm{2}},\mathit{\varphi })=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}})}}\\ {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{w_{\mathrm{1}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}}}+{\displaystyle \frac{w_{\mathrm{2}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}\mathit{\rho }{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{U}}}_{\mathrm{2}}\cos ({\mathit{\varphi }}_{\mathrm{1}}-{\mathit{\varphi }}_{\mathrm{2}}+\mathit{\psi })}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}}\right]\\ {\displaystyle }& {\displaystyle }+\ln I_{\mathrm{0}}\left(\sqrt{A^{\mathrm{2}}+B^{\mathrm{2}}+\mathrm{2}AB\cos \mathit{\lambda }}\right)\\ \text{(A2)}& {\displaystyle }& {\displaystyle }+C\cos \mathit{\nu }\cos \mathit{\lambda }+\ln \cosh \left(C\sin \mathit{\nu }\sin \mathit{\lambda }\right)\mathrm{d}\mathit{\lambda },\end{array}$$

with

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }A={\displaystyle \frac{w_{\mathrm{1}}\sqrt{a_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{1}}^{\mathrm{2}}}}{{\mathit{\sigma }}_{\mathrm{1}}}},\text{(A4)}& {\displaystyle }& {\displaystyle }B={\displaystyle \frac{w_{\mathrm{2}}\sqrt{a_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{2}}^{\mathrm{2}}}}{{\mathit{\sigma }}_{\mathrm{2}}}},\text{(A5)}& {\displaystyle }& {\displaystyle }C={\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}{\displaystyle \frac{w_{\mathrm{1}}{w}_{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}}},\end{array}$$

and the integral is evaluated numerically. Equation (A1) is obtained from Eq. (26) using the fact that

$$\begin{array}{}\text{(A6)}& I_k\left(x\right)I_k\left(y\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{I}_{\mathrm{0}}\left(\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}+\mathrm{2}xy\cos \mathit{\varphi }}\right)\cos k\mathit{\varphi }\mathrm{d}\mathit{\varphi }\end{array}$$

(which follows from Eq. 33) and use of Eq. (24). When w_i∕σ_i becomes large, numerical evaluations of the Bessel function in Eq. (A1) become unreliable. For the present computations using Matlab, values of Inf occur in such cases. This problem was not solved by using the approximation Eq. (10) when the argument of the Bessel function is large, as this approximation is not sufficiently accurate.

Appendix B: Bivariate goodness-of-fit test

Goodness-of-fit of the bivariate distributions considered was assessed as follows. For the speed dataset w_j,n, $j=\mathrm{1},\mathrm{2}$ and $n=\mathrm{1},\mathrm{\dots },N$, evenly spaced quantiles $q_{j,i}=i/M$, $i=\mathrm{0},\mathrm{\dots },M$ for the marginals are estimated. The quantiles $q_{j,\mathrm{0}}=\mathrm{0}$ and $q_{j,M}=\mathrm{1}$ are estimated respectively as 0.9 times the smallest observed value and 1.1 times the largest observed value. The number of pairs of observations falling simultaneously into all pairs of quantiles are computed:

$$\begin{array}{ll}\text{(B1)}& {\displaystyle }& {\displaystyle }{f}_{kl}={\displaystyle }& {\displaystyle }\sum _{n=\mathrm{1}}^N\mathbf{\text{1}}\left[\left(w_{\mathrm{1},n}\in (q_{\mathrm{1},k},q_{\mathrm{1},(k+\mathrm{1})}]\right)\cap \left(w_{\mathrm{2},n}\in (q_{\mathrm{2},l},q_{\mathrm{2},(l+\mathrm{1})}]\right)\right],\end{array}$$

where 1(⋅) is the indicator function. The pdf with maximum likelihood parameters θ is then integrated to obtain the expected number of observations in these intervals

$$\begin{array}{}\text{(B2)}& g_{kl}=N\underset{q_{\mathrm{1},k}}{\overset{q_{\mathrm{1},(k+\mathrm{1})}}{\int }}\underset{q_{\mathrm{2},l}}{\overset{q_{\mathrm{2},(l+\mathrm{1})}}{\int }}p(w_{\mathrm{1}},w_{\mathrm{2}};\mathit{\theta })\mathrm{d}{w}_{\mathrm{1}}\mathrm{d}{w}_{\mathrm{2}}\end{array}$$

and the test statistic A is computed:

$$\begin{array}{}\text{(B3)}& A={\displaystyle \frac{\mathrm{1}}{M^{\mathrm{2}}}}\sum _{k,l=\mathrm{1}}^M|f_{kl}-g_{kl}|.\end{array}$$

Any elements g_kl which take the value of Inf (because of numerical difficulties in evaluating the Bessel function for large arguments; cf. Appendix A) are excluded from the calculation of the test statistic A.

After computation of A from the observations, an ensemble of B random realizations of length N from $p(w_{\mathrm{1}},w_{\mathrm{2}};\mathit{\theta })$ is generated and the corresponding ${\stackrel{\mathrm{̃}}{A}}_k$, $k=\mathrm{1},\mathrm{\dots },B$ values of the test statistic are generated in the manner described above. The p-value of the null hypothesis that the observations are drawn from the specified distribution is finally computed as the fraction of ${\stackrel{\mathrm{̃}}{A}}_k$ values falling above A. Throughout this analysis, we use M=20 and B=250 (unless otherwise noted). This goodness-of-fit test assumes independence of the random draws $(w_{\mathrm{1}n},w_{\mathrm{2}n}),(w_{\mathrm{1}m},w_{\mathrm{2}m})$, m≠n. To minimize the effect of serial dependence in data, in this study we subsample the datasets considered with a sampling interval sufficiently large to balance reducing serial dependence with maintaining sample size.

A second goodness-of-fit test proposed by McAssey (2013) was also considered, but the statistical power was found to be lower than the test described above for the distributions considered in this study.