Wind energy applications including wind resource assessment, wind power forecasting, and wind plant optimization require high-resolution mesoscale simulations. High resolution mesoscale simulations are essential for accurate characterization of atmospheric flows over heterogeneous land use and complex terrain. Under such conditions, the assumption of grid-cell homogeneity, used in one-dimensional planetary boundary layer (1D PBL) parameterizations, breaks down. However, in most numerical weather prediction (NWP) models, boundary layer turbulence is parameterized using 1D PBL parameterizations. We have therefore developed a three-dimensional (3D) PBL parameterization to better account for horizontal flow heterogeneities. We have implemented and tested the 3D PBL parameterization in the Weather Research and Forecasting (WRF) numerical weather prediction model. The new parameterization is validated using observations from the Wind Forecast Improvement 2 (WFIP 2) project and compared to 1D PBL results.