The goal of this study was to develop a statistical bivariate wind speed-wind shear model (WSWS). The development of WSWS is based on near surface wind speed data available from 397 measurement stations distributed over Germany, as well as on ERA-Interim reanalysis wind speed data available in 1000 m above ground level (a.g.l.). These data were used (1) to calculate empirical distributions of wind speed in 1000 m a.g.l., (2) empirical distributions of the wind shear exponent, and (3) to fit theoretical distributions to the empirical wind speed and wind shear exponent distributions. It was found that the four parameter Johnson SB distribution reproduces the shape of the wind speed in 1000 m a.g.l. empirical distributions best. The four parameter Dagum distribution provided good fits to the empirical wind shear distributions. The parameterized wind speed and wind shear marginal distributions were then linked by 16 joint copulas. Goodness-of-fit evaluation of the joint copulas demonstrates that the Gaussian-Gaussian copula reproduces the empirical bivariate wind speed-wind shear distribution most accurately. By using WSWS it is possible to continuously calculate the wind speed probability density function in hub heights between 10 m a.g.l. and 200 m a.g.l. This allows WSWS to be applied to virtually any power curve for computing the wind energy yield and capacity factor in the analyzed hub height range. A one-time site-specific parametrization of WSWS is sufficient for a comprehensive height-dependent exploitation of the available wind resource. (C) 2017 Elsevier Ltd. All rights reserved.