Optimizing the Weather Research and Forecasting ( WRF ) Model for Mapping the Near-Surface Wind Resources over the Southernmost Caribbean Islands of Trinidad and Tobago

Numerical wind mapping is currently the wind power industry’s standard for preliminary assessments for sites of good wind resources. Of the various available numerical models, numerical weather prediction (NWP) models are best suited for modeling mesoscale wind flows across small islands. In this study, the Weather Research and Forecast (WRF) NWP model was optimized for simulating the wind resources of the Caribbean islands of Trinidad and Tobago in terms of spin-up period for developing mesoscale features, the input initial and boundary conditions, and the planetary boundary layer (PBL) parameterizations. Hourly model simulations of wind speed and wind direction for a one-month period were compared with corresponding 10 m level wind observations. The National Center for Environmental Prediction (NCEP)-Department of Energy (DOE) reanalysis of 1.875◦ horizontal resolution was found to be better suited to provide initial and boundary conditions (ICBCs) over the 1◦ resolution NCEP final analysis (FNL); 86% of modeled wind speeds were within ±2 m/s of measured values at two locations when the coarse resolution NCEP reanalysis was used as compared with 55–64% of modeled wind speeds derived from FNL. Among seven PBL schemes tested, the Yonsei University PBL scheme with topographic drag enabled minimizes the spatial error in wind speed (mean bias error +0.16 m/s, root-mean-square error 1.53 m/s and mean absolute error 1.21 m/s) and is capable of modeling the bimodal wind speed histogram. These sensitivity tests provide a suitable configuration for the WRF model for mapping the wind resources over Trinidad and Tobago, which is a factor in developing a wind energy sector in these islands.


Introduction
The economy and gross domestic product (GDP) of Trinidad and Tobago, the southernmost islands of the Caribbean, are hinged on the energy sector, with as much as 80% of export revenues coming from the petrochemical sector [1].The indigenous natural gas supplies are used as fuel for electricity on the islands, enabling consumers to pay the second lowest residential electricity tariff in the Caribbean region (US $0.06/kWh) [2].In terms of availability of resources, solar and wind energies are all suitable forms of renewable energy for electricity generation on the islands and will enable a diversification of the energy mix.Wind energy is expected to play a major role for bulk Energies 2017, 10, 931; doi:10.3390/en10070931www.mdpi.com/journal/energieselectricity generation in the twin-island state of Trinidad and Tobago [3] and has been shown to be more cost-competitive than solar energy in other regions [4].
The very first step toward wind farm development is to rank potential sites based on wind resources, land availability, location and stability of the electricity grid, environmental impacts, and other factors.This step therefore requires a good understanding of the spatial wind resources and availability.To aid in renewable energy decisions, the International Renewable Energy Agency (IRENA) has made available wind maps for all countries at a horizontal resolution of 5 km × 5 km from Vaisala [5].In addition, 250 m horizontal resolution wind maps [5] constructed by Denmark Technical University (DTU) using the Wind Atlas Analysis and Application program have also been provided by IRENA.Vaisala's maps were generated using a prognostic numerical weather prediction (NWP) model [6] and DTU's with a diagnostic wind flow adjustment model [7].Of these two types of numerical models, NWPs are better suited for modelling local and regional wind circulation flows in complex coastal topography, including the prevalent sea and land breezes around islands.Furthermore, NWPs have been used for preliminary wind resource assessment in several regions globally [8][9][10][11][12][13][14][15] with few studies for the Caribbean region.
Several studies [15][16][17][18] indicate that NWP model configuration for simulating wind flow is region dependent.However, only one specific model configuration of a NWP was used in producing Vaisala's global wind maps.NWP studies for the Caribbean region have focused on a NWP's ability and sensitivity in predicting short-term heavy rainfall events in Puerto Rico, Barbados and Guyana [19][20][21] and hurricane evolution and associated storm surge [22][23][24].Therefore, this study intends to optimize a NWP model output to produce more accurate model winds for one of the Caribbean small island states.In wind resource assessments, NWP model winds are used to determine mean wind speed maps, wind speed distributions, and wind power density which vary with the cube of the wind speed.As relative errors in wind power density are three times the relative errors in model wind speed, accurate predictions of wind speed will result in reliable maps of wind power density.
The accurate simulation of low-level winds at high resolutions in the NWP context is dependent on a few key components of NWP models.NWP models simulate the atmosphere as a moving fluid as well as the chemical reactions within the atmosphere on user defined grid scales.They also account for the influence of sub-grid physical and chemical processes on grid-resolved motions through several types of parameterizations.The components of the NWP model that influence the accuracy of the simulations include the initial and boundary atmospheric conditions (reanalysis data) (e.g., [25]), the spin-up time the model takes to generate mesoscale information that is absent from the coarse boundary conditions (e.g., [26]), and the parameterization of the influence of turbulence on wind velocity components (e.g., [18,27]).
Several institutions produce reanalysis data across the globe by assimilating measured data from available sources into general circulation models (GCMs).The result of this is a spatially and temporally homogeneous gridded data set of atmospheric conditions [28].Among the reanalysis products that are available, the most popular ones that have been successfully downscaled for numerical wind mapping and wind forecasting include the National Center for Environmental Protection/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (e.g., [29,30]), and the European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis (ERA-Interim) (e.g., [31,32]).Few studies have evaluated reanalysis products for the Caribbean region although the National Center for Environmental Protection/Department of Energy (NCEP/DOE) reanalysis, an updated version of the NCEP/NCAR reanalysis, has been shown to provide realistic representations of the large-scale wind power in the Caribbean region [33,34].In addition to traditional reanalysis, previous studies on NWP wind resource mapping have commonly used the NCEP final analysis (FNL) [10,35] which differs from the traditional reanalysis in that reanalysis assimilates a higher quantity of observed data as they are generated with non-real time observational data sets [25].The FNL analysis is subjected to changes in the operational configuration of the global model [36].
Energies 2017, 10, 931 3 of 23 Large-scale atmospheric conditions do not contain realistic vertical motions associated with mesoscale coastal circulations.At initial time, mesoscale information are non-existent.NWP models add smaller-scale mesoscale information than that provided by the input large-scale (synoptic) atmospheric conditions and can do so via two methods: data assimilation and spin-up periods.The assimilation of measurements into the NWP for regions of sparse data sets could degrade the performance of the model rather than improve it because the mesoscale and local scale may not be adequately represented due to the large spatial extent between measurement stations.Alternatively, the model could spin-up these features during the initial hours of the forecasts/hindcasts via a spin-up period to develop small-scale information on scales less than the resolution of the initial conditions [37].The spin-up period is therefore a period during which the model approaches physical and dynamical equilibrium.This spin-up period is discarded and not used in subsequent analyses.Typical spin up times vary from 2 to 3 [38,39], 6 [16,40], 12 h [41,42] to 24 h [31] and could extend to days [37].
Of the various physical and chemical parameterizations, the coupled planetary boundary layer (PBL) and surface layer (SL) processes have been found to consistently influence boundary layer winds the most [43][44][45].As the PBL-SL parameterization account for the influence of turbulence on wind velocity components, the selection of an appropriate PBL-SL (hereafter PBL) scheme is required to better estimate mean wind speeds in the PBL since small relative errors in the wind speed are magnified three times in the relative errors in wind power densities.There is a need to determine the PBL scheme that performs best for each region.Different PBL schemes provided an overall best performance for Denmark [16], the North Sea [17], Southern Spain [43] and Texas [46].
Therefore, a set of NWP model configurations has to be determined for mapping the spatial wind resources.In this study, we configure the WRF NWP for mapping the wind resources of the island state of Trinidad and Tobago, the southernmost islands in the Caribbean, as this model has been tested for mapping the wind resources in several countries [17,25,[47][48][49].The objectives of this study were to determine the spin-up period for developing mesoscale features as well as assess which one of two input atmospheric data sets and seven PBL parameterizations available to the WRF model would optimize the accuracy of simulating the wind resources for the area of interest.In Section 2 we provide descriptions of the area of interest, the initial configuration for the WRF model, the sensitivity tests to input atmospheric data sets as well as those to PBL schemes, and the in-situ data that were used to compute error metrics and in determining model configuration.The results of the sensitivity tests are described and discussed in Section 3 and conclusions are made in Section 4.

Area of Interest
The numerical domain on which the optimization of the WRF NWP model was performed is shown in Figure 1.It consists of two telescoping nests centered on 10.69 • N and 61.47 • W with horizontal resolutions of 25 km on the outermost domain (d01), 5 km on the first nest (d02), and 1 km on the innermost nest (d03).d01 consists of 91 × 61 grid-points and covers the Eastern Caribbean islands, part of South America and portions of the Caribbean Sea and the Atlantic Ocean.The first nest d02 comprises 151 × 111 grid-points and contains Trinidad and Tobago and the north-eastern part of Venezuela.The second nest d03 encompasses Trinidad and Tobago only, the area of interest, and is spanned by 331 × 196 grid-points.The numerical grid was configured with 34 vertical levels which follow the terrain.The first six levels of the grid above ground were approximately 14 m, 42 m, 76 m, 117 m, 191 m, and 277 m which are close to contemporary wind turbine hub heights at approximately 50 m, 60 m, 80 m, 100 m and 120 m.The vertical grid levels through a latitudinal cross-section (10.72 • N) are shown in Figure 2a.The highest point in Trinidad and Tobago in nest d03 is located along this latitude.Figure 2b shows the first nine levels of the grid with the first vertical grid level indicated in blue and the terrain profile in red.The heights shown in Figure 2a,b are the model heights above sea level.The heights of these vertical levels may vary slightly from one location to another.Figure 2c      Energies 2017, 10, 931 5 of 23

Overview of Method and Data
The large-scale atmospheric data sets used in this study are described in Section 2.2.1.In Section 2.2.2, the in-situ data and error metrics against which the high-resolution simulations were tested to determine the suitability of the input atmospheric conditions are presented.The method by which the spin-up period is determined is outlined in Section 2.2.3.In Section 2.2.4, the initial configuration for the WRF model is specified and the PBL-SL schemes that were tested to determine a suitable one for future wind mapping studies are described briefly in Section 2.2.5.

Reanalysis Data for NWP Model Initialization
The WRF model ingests, as input, large-scale atmospheric conditions that were used to initialize the models prior to time integration.These conditions are modified by the WRF model for local topographic conditions.Two atmospheric data sets were tested as initial and boundary conditions.As the 10 m winds of the National Center for Environmental Prediction/Department of Energy (NCEP/DOE) reanalysis (1.875 • × 1.875 • resolution) were found to provide reasonable estimates of mean wind speeds and wind power densities for the Caribbean [32], we tested this reanalysis as well as the higher resolution NCEP's global tropospheric final analysis (FNL) (1 • × 1 • ) for a randomly selected one-month period, April 2003.Hourly simulations with the WRF model for one month should be able to capture the diurnal changes in the wind.Both NCEP-DOE reanalysis (hereafter NCEP reanalysis) and FNL are available every six hours (00, 06, 12 and 18 UTC (coordinated universal time)).The NCEP-DOE reanalysis is based on the climate data assimilation system (CDAS2) which is an older version of the global forecast system (GFS) [28].The FNL analysis is a product of the global data assimilation system (GDAS2) for which the GFS is the general circulation model [28].WRF is reinitialized from a "cold start" at 1200 UTC daily and run for 36 h for spin-up tests in the first instance for each day of April 2003.

In-Situ Data and Error Metrics
The WRF modeled 10 m wind speed and wind direction were compared to hourly in-situ observations at the two anemometric stations known to adhere to observation standards of the World Meteorological Organization (WMO): Piarco, Trinidad (10.617 • N, 61.350 • W) and Crown Point, Tobago (11.150 • N, 60.833 • W).Both stations have low elevations, 12 m and 15 m at Crown Point and Piarco, respectively.The Meteorological Service of Trinidad and Tobago is responsible for climatic data measurement and collection and these data are available from the US National Climatic Data Center [50].The predicted 10 m wind speeds at Crown Point and Piarco were determined by bilinearly interpolating the meridional (u) and zonal (v) wind components of the four surrounding grid-points to each station's latitude and longitude.
For N predicted and observed scalar pairs at forecast hour i, (U P,i , U O,i ), the following statistics were calculated [27,[51][52][53]: mean predicted wind speed U P , mean observed wind speed U O , mean absolute error (MABE), root-mean-square error (RMSE), mean bias error (MBE) and the standard error in the wind speed (STDE).In addition, the Pearson correlation coefficient (R 2 ) and the index of agreement (IOA) between the predicted hourly wind speeds and their corresponding observed values were computed.These error statistics are described in Appendix A.

Spin-Up Period
The input analyzed data contain kinetic energy at the synoptic scales and not at the mesoscale.Spin-up periods are therefore determined from the generation of energy over the mesoscale by computing one dimensional kinetic energy spectra of the zonal and meridional wind speed components on the innermost grid d03 at several hours into the runs [26].Wind components were detrended to make the data series periodic prior to computing the normalized one dimensional kinetic energy spectra along the west-east horizontal grid lines neglecting 15 grid-points from the west and east Energies 2017, 10, 931 6 of 23 boundaries of the innermost domain [26].Spectra were also averaged over the north-south extent again neglecting 15 grid-points from the south and north boundaries.This procedure was repeated for 0, 6, 12, 18 and 24 h into the 36-h forecasts.

Initial Configuration for the Weather Research and Forecasting (WRF) Model
The physical parameterizations of the Advanced Research Weather Research and Forecasting (WRF) model, version 3.5.1 [54], that were used for the first experiment were the rapid radiative transfer model (RRTM) for longwave radiation, the Dudhia scheme for shortwave radiation, the Yonsei University planetary boundary layer (PBL) scheme which works in conjunction with the Monin-Obukhov surface-layer (SL) parameterization, the thermal diffusion land surface model, and the WRF single-moment 3 (WSM 3) scheme for microphysics.The Betts-Miller-Janjic (BMJ) scheme for cumulus parameterization was applied to the outermost (25 km resolution) domain only, while cloud formation was explicitly resolved on the nests (5 km and 1 km resolutions).Rama Rao et al. [21] had found that the combination of the BMJ cumulus and the WSM 3 microphysics schemes were suitable for reproducing intense rainfall events over Guyana (~6 • N).Rama Rao's simulation domain had included the Caribbean's southernmost island of Trinidad [21].The WRF land cover characterization (LCC) at 30 arc-second (~900 m) resolution provided land use globally in 24 categories.The WRF LCC is a version of the United States geological survey/global land cover characterization (USGS/GLCC) and works in conjunction with the USGS 30 arc-second elevation data (GTOPO30) by WRF's preprocessor (WPS) to generate input topographical data for WRF.Atmospheric boundary conditions at the outermost grid were updated every 6 h without observational or grid nudging as the integration period is short enough to prevent excessive model drift from reality.In addition to atmospheric fields, skin temperatures in the analyzed/reanalysis were used to provide sea surface temperatures (SST).SST once initialized each day was held constant throughout the integration period.WRF was run on time steps of 120 s on the parent domain.

Planetary Boundary Layer (PBL) Schemes
The WRF model has several PBL schemes, each PBL scheme having different ways of representing the vertical mixing of turbulence in the PBL.We vary the PBL scheme to determine which one best represents the mean surface wind speeds.The Yonsei University (YSU) scheme [55] was used in the first experiment.This PBL scheme is a first-order, non-local closure scheme.The first-order scheme expresses the turbulent fluxes in terms of an eddy diffusivity and the vertical gradient of the mean variable and includes a counter-gradient term.In the first simulation experiment, the topographic correction for extra drag from subgrid topography in the YSU scheme is not enabled.This experiment is referred to as YSU-NO-TOPO.The second experiment used YSU PBL but with topographic correction enabled (YSU-TOPO).The second PBL scheme tested was the asymmetric convective model (ACM2) [56].The ACM2 is also a first-order closure with a combination of local and non-local mixing approach and accounts for non-local mixing between adjacent layers.The Mellor-Yamada-Janjic (MYJ) [57], Mellor-Yamada-Nakanishi-Niino (MYNN 2.5) [58,59], Quasi-normal scale elimination (QNSE) [60], and the University of Washington scheme (UW-TK) [61]) are classified as turbulent kinetic energy (TKE) closure schemes that allow only local transport.TKE closure schemes include a prognostic equation for the TKE unlike the first-order scheme.The total energy mass-flux (TEMF) [62] scheme is also a TKE closure scheme that also incorporates a mass flux for non-local vertical mixing.The eddy diffusivity in TKE closure schemes are functions of the mixing length for which different formulations exist depending on the parameterization.Most TKE schemes in WRF with the exception of the MYNN 2.5 are 1.5 order schemes; the MYNN 2.5 is a 2.5 order closure scheme.

Spin-Up Period
Kinetic energy spectra were determined for each day in April 2003 at every 6 h interval of the integration period.Figures 3 and 4 show the corresponding kinetic energy spectra for only two days, 2 April 2003 and 6 April 2003, as other days show similar spectra.Figure 3 spectra were derived from model simulations that ingested FNL data as initial while those in Figure 4 were from simulations initialized with NCEP reanalysis.The spin-up period is determined by comparing the slope of the kinetic energy spectra with the κ −5/3 (κ is the wavenumber) dependence which has been observed at the mesoscale scales [63][64][65].At initialization, that is, zero hours into the forecast, no smaller-scale information is present and the energy content is low (Figures 3 and 4).Substantial mesoscale information develops by 6 h into the simulation.Spectra for 12 h and 18 h are similar to the 6 h spectra for 2 April 2003.However, for 6 April 2003 the 6 h forecast has not developed sufficiently; by 12 h into the run mesoscale information seems to have fully developed.We therefore chose the first 12 h as the spin-up period for generation of mesoscale information.Thus, the simulation of one diurnal period is a 36 h run with the first 12 h discarded for the development of mesoscale information.
The WRF model develops sufficient mesoscale information 12 h from model initialization thereby providing a spin-up time configuration of 12 h.In addition, the spin-up period of 12 h is in keeping with the corresponding timescales associated with typical mesoscale features.We note that the modelling of some days require a shorter spin-up time of 6 h, and that spin-up time may be dependent on the large-scale atmospheric conditions driving mesoscale phenomena.Therefore, errors between the simulated wind speeds and corresponding observations may be introduced by using a fixed spin-up period of 12 h for all synoptic conditions.However, such errors may be smoothed out when deriving average wind statistics relevant for wind maps.We also found that Skamarock's spectral method [22] in which the kinetic energy spectra of model simulations were compared with the typical spectra of mesoscale features is useful for model configuration for regions with limited data sets.
2 April 2003 and 6 April 2003, as other days show similar spectra.Figure 3 spectra were derived from model simulations that ingested FNL data as initial while those in Figure 4 were from simulations initialized with NCEP reanalysis.The spin-up period is determined by comparing the slope of the kinetic energy spectra with the / ( is the wavenumber) dependence which has been observed at the mesoscale scales [63][64][65].At initialization, that is, zero hours into the forecast, no smaller-scale information is present and the energy content is low (Figures 3 and 4).Substantial mesoscale information develops by 6 h into the simulation.Spectra for 12 h and 18 h are similar to the 6 h spectra for 2 April 2003.However, for 6 April 2003 the 6 h forecast has not developed sufficiently; by 12 h into the run mesoscale information seems to have fully developed.We therefore chose the first 12 h as the spin-up period for generation of mesoscale information.Thus, the simulation of one diurnal period is a 36 h run with the first 12 h discarded for the development of mesoscale information.
The WRF model develops sufficient mesoscale information 12 h from model initialization thereby providing a spin-up time configuration of 12 h.In addition, the spin-up period of 12 h is in keeping with the corresponding timescales associated with typical mesoscale features.We note that the modelling of some days require a shorter spin-up time of 6 h, and that spin-up time may be dependent on the large-scale atmospheric conditions driving mesoscale phenomena.Therefore, errors between the simulated wind speeds and corresponding observations may be introduced by using a fixed spin-up period of 12 h for all synoptic conditions.However, such errors may be smoothed out when deriving average wind statistics relevant for wind maps.We also found that Skamarock's spectral method [22] in which the kinetic energy spectra of model simulations were compared with the typical spectra of mesoscale features is useful for model configuration for regions with limited data sets.

General Observations for Hourly Simulations of Wind Speed
The modeled hourly time series for the first fourteen days of April 2003 are shown in Figures 5 and 6 for Crown Point and Piarco, respectively.The WRF model's configuration with the YSU PBL and with both NCEP and FNL initial conditions (ICs) is capable of reproducing the diurnal behavior of the wind speed at both stations.However, the configurations seem to underestimate the diurnal maxima at Crown Point.The model overestimates the nighttime low wind speeds at both stations.The model underestimates wind speeds that exceed 8 m/s at Piarco.In general, the model configuration and the input atmospheric data capture the timing of the daily maxima and minima moreso at Piarco than Crown Point.There are however, some timing problems at Crown Point on 4 April (~72 h) and 12 April (~280 h).

General Observations for Hourly Simulations of Wind Speed
The modeled hourly time series for the first fourteen days of April 2003 are shown in Figures 5 and 6 for Crown Point and Piarco, respectively.The WRF model's configuration with the YSU PBL and with both NCEP and FNL initial conditions (ICs) is capable of reproducing the diurnal behavior of the wind speed at both stations.However, the configurations seem to underestimate the diurnal maxima at Crown Point.The model overestimates the nighttime low wind speeds at both stations.The model underestimates wind speeds that exceed 8 m/s at Piarco.In general, the model configuration and the input atmospheric data capture the timing of the daily maxima and minima moreso at Piarco than Crown Point.There are however, some timing problems at Crown Point on 4 April (~72 h) and 12 April (~280 h).The diurnal behavior at Crown Point may be more difficult for WRF to capture because it is a coastal station with sea-land breeze circulations while Piarco is an inland station.However it is more likely that the diurnal variation at Piarco was better modeled than Crown Point because the data values for this station may have been assimilated during the production of the analysis/reanalysis data.

Error Statistics for FNL vs NCEP as Initial and Boundary Conditions
The diurnal variability of the error statistics are shown in Figures 7 and 8.At Crown Point, the hourly averaged errors for model simulations using NCEP ICs are either smaller or similar to the model errors obtained using FNL ICs.In addition, error statistics at Crown Point are smaller for 12 UTC to 20 UTC (about 7 A.M. to 3 P.M. LST (local standard time)) than other times of the day.Daytime convective conditions may be better reproduced at Crown Point than nighttime stable conditions.However, the observed in-situ wind data were taken at only one height and thus, no verifications can be made for stability conditions or wind shear.At Piarco, model errors obtained on using NCEP as atmospheric input are less than those from FNL for about 12 h.
The average error statistics for the two stations for April 2003 are shown in Table 1.Error statistics (MBE, MABE, RMSE) that were computed using WRF model output at the 1 km resolution and station data indicate that using NCEP and FNL produce similar wind speed values with NCEP possibly producing a greater number of more realistic hourly predictions of the mean wind speed.The bias errors with NCEP are a small negative value for Piarco and a positive at Crown Point.The mean error statistics of MABE and RMSE are similar for NCEP and FNL input conditions.The relatively low R 2 values at Piarco and Crown Point indicate that weak linear relationship exist between the simulated winds and the observations.As the R 2 statistic is not sensitive to additive and proportional differences between simulations and observations [52], we also computed the IOA.The relatively high IOA values (0.73 at Crown Point and 0.88 at Piarco) indicate that the phase relationships between the simulated and in-situ data are reasonably well reproduced with fewer timing errors at Piarco than Crown Point.The diurnal behavior at Crown Point may be more difficult for WRF to capture because it is a coastal station with sea-land breeze circulations while Piarco is an inland station.However it is more likely that the diurnal variation at Piarco was better modeled than Crown Point because the data values for this station may have been assimilated during the production of the analysis/reanalysis data.

Error Statistics for FNL vs NCEP as Initial and Boundary Conditions
The diurnal variability of the error statistics are shown in Figures 7 and 8.At Crown Point, the hourly averaged errors for model simulations using NCEP ICs are either smaller or similar to the model errors obtained using FNL ICs.In addition, error statistics at Crown Point are smaller for 12 UTC to 20 UTC (about 7 A.M. to 3 P.M. LST (local standard time)) than other times of the day.Daytime convective conditions may be better reproduced at Crown Point than nighttime stable conditions.However, the observed in-situ wind data were taken at only one height and thus, no verifications can be made for stability conditions or wind shear.At Piarco, model errors obtained on using NCEP as atmospheric input are less than those from FNL for about 12 h.The average error statistics for the two stations for April 2003 are shown in Table 1.Error statistics (MBE, MABE, RMSE) that were computed using WRF model output at the 1 km resolution and station data indicate that using NCEP and FNL produce similar wind speed values with NCEP possibly producing a greater number of more realistic hourly predictions of the mean wind speed.The bias errors with NCEP are a small negative value for Piarco and a positive at Crown Point.The mean error statistics of MABE and RMSE are similar for NCEP and FNL input conditions.The relatively low R 2 values at Piarco and Crown Point indicate that weak linear relationship exist between the simulated winds and the observations.As the R 2 statistic is not sensitive to additive and proportional differences between simulations and observations [52], we also computed the IOA.The relatively high IOA values (0.73 at Crown Point and 0.88 at Piarco) indicate that the phase relationships between the simulated and in-situ data are reasonably well reproduced with fewer timing errors at Piarco than Crown Point.
Table 1.Model and error statistics of WRF's 10 m level wind speed and wind direction at Piarco and Crown Point using FNL and NCEP-DOE analyses as initial/boundary conditions.WRF was run with two nests and modeled wind speeds were found on the inner-most domain which has horizontal resolution of 1 km (DO3 on Figure 1).Error statistics with the "dir" subscript refers to the error in wind direction.Error statistics were computed for 720 h of simulated winds.While wind directional RMSE statistics (RMSE dir = 17-27 • ) are small compared with those found for other regions such as Portugal (37-43 • ) [21], the bias in wind direction are large, 9-15 • in this study and −2-0.5 • in Carvalho et al. [21] for Portugal.The mean absolute error and mean bias error in wind directions (MAE dir , MBE dir ) are similar for both sites on using the two input atmospheric data sets.However, NCEP produces a smaller RMSE error in the wind direction (RMSE dir ) at Piarco.The small wind directional errors may be due to the low dispersive nature of the actual wind directional distributions at both Crown Point and Piarco (Figure 9, red bars) which are reproduced by the WRF model although the bimodality of the wind direction distribution at Piarco (Figure 9b) is not.The WRF-NCEP and WRF-FNL also reproduced the greater percentage of east-south-easterly wind at Crown Point (Figure 9a) than Piarco (Figure 9b).

Statistic
While wind directional RMSE statistics ( = 17-27°) are small compared with those found for other regions such as Portugal (37-43°) [21], the bias in wind direction are large, 9-15° in this study and −2-0.5° in Carvalho et al. [21] for Portugal.The mean absolute error and mean bias error in wind directions ( , ) are similar for both sites on using the two input atmospheric data sets.However, NCEP produces a smaller RMSE error in the wind direction ( ) at Piarco.The small wind directional errors may be due to the low dispersive nature of the actual wind directional distributions at both Crown Point and Piarco (Figure 9, red bars) which are reproduced by the WRF model although the bimodality of the wind direction distribution at Piarco (Figure 9b) is not.The WRF-NCEP and WRF-FNL also reproduced the greater percentage of east-south-easterly wind at Crown Point (Figure 9a) than Piarco (Figure 9b).Comparisons of wind speed histograms of simulated wind speed with in-situ data (Figures 10 and 11) show that the frequency of wind speeds in the 2-3 m/s and the 3-4 m/s wind speed intervals of the in-situ data are over-predicted at both Piarco and Crown Point.At Piarco Comparisons of wind speed histograms of simulated wind speed with in-situ data (Figures 10  and 11) show that the frequency of wind speeds in the 2-3 m/s and the 3-4 m/s wind speed intervals of the in-situ data are over-predicted at both Piarco and Crown Point.At Piarco (Figure 11), the frequency of wind speeds greater than 6 m/s is reasonably well predicted while at Crown Point (Figure 10) this frequency is under-predicted.We do not compute chi-squared statistics to compare modeled distributions with that of the in-situ data since these may not be independent data sets as some of the station data may have been assimilated during the production of the reanalysis and analysis coarse resolution initial and boundary conditions.The wind speed frequency distribution at the near-surface inland station at Piarco is better predicted than that at the coastal station at Crown Point.Since the frequency of wind speeds >6 m/s at these two stations are under-predicted, the estimates of the WRF model using either FNL or NCEP atmospheric conditions would most likely produce conservative estimates of the wind resource.
In lieu of a direct comparison between the modeled and in-situ wind speed distributions, the wind speed error distributions were compared (Figures 12 and 13).The error distribution adds information that is not accounted for in the mean bias error.The mean bias error with NCEP is only slightly better than FNL, 1.13 m/s versus 1.32 m/s at Crown Point, and −0.01 m/s and 0.05 m/s at Piarco.
However, the wind speed error distributions differentiate between large hourly negative and positive errors which are averaged to determine the MBE statistic.The error distribution/histogram determines how often the NWP make unacceptable errors.Model results at Crown Point with NCEP as input data produced 86% of modeled hourly values within ±2 m/s of the actual values while those with FNL input data produced 64% (Figure 12).The percentage of hourly modeled wind speeds within ±1 m/s of the actual values were 55% with NCEP input data and 34% for FNL input data.At Piarco, WRF results using either FNL or NCEP as input have approximately 87% of modeled hourly values within ±2 m/s of the actual wind speed and 56% within ±1 m/s (Figure 13).
(Figure 11), the frequency of wind speeds greater than 6 m/s is reasonably well predicted while at Crown Point (Figure 10) this frequency is under-predicted.We do not compute chi-squared statistics to compare modeled distributions with that of the in-situ data since these may not be independent data sets as some of the station data may have been assimilated during the production of the reanalysis and analysis coarse resolution initial and boundary conditions.The wind speed frequency distribution at the near-surface inland station at Piarco is better predicted than that at the coastal station at Crown Point.Since the frequency of wind speeds >6 m/s at these two stations are under-predicted, the estimates of the WRF model using either FNL or NCEP atmospheric conditions would most likely produce conservative estimates of the wind resource.In lieu of a direct comparison between the modeled and in-situ wind speed distributions, the wind speed error distributions were compared (Figures 12 and 13).The error distribution adds information that is not accounted for in the mean bias error.The mean bias error with NCEP is only Crown Point (Figure 10) this frequency is under-predicted.We do not compute chi-squared statistics to compare modeled distributions with that of the in-situ data since these may not be independent data sets as some of the station data may have been assimilated during the production of the reanalysis and analysis coarse resolution initial and boundary conditions.The wind speed frequency distribution at the near-surface inland station at Piarco is better predicted than that at the coastal station at Crown Point.Since the frequency of wind speeds >6 m/s at these two stations are under-predicted, the estimates of the WRF model using either FNL or NCEP atmospheric conditions would most likely produce conservative estimates of the wind resource.In lieu of a direct comparison between the modeled and in-situ wind speed distributions, the wind speed error distributions were compared (Figures 12 and 13).The error distribution adds information that is not accounted for in the mean bias error.The mean bias error with NCEP is only In summary, the coarser NCEP-DOE reanalysis (1.875 • horizontal resolution) provides more accurate estimates of the wind speeds than the FNL (1 • resolution) analysis in terms of minimizing mean wind speed errors and producing more symmetric wind speed error distributions for the two measurement stations in Trinidad and Tobago.Sensitivity tests for input and boundary conditions found that the coarser NCEP-DOE reanalysis produced more accurate predictions of mean wind speed, smaller wind speed biases, and a greater percentage of wind speed errors within acceptable limits than the higher resolution FNL analysis.One possible explanation for NCEP-DOE providing more accurate estimates of the mean wind speeds over the FNL may be the differences in the data assimilation schemes that were used to produce the analyses.In future, as more wind data sets in Energies 2017, 10, 931 13 of 23 tropical locations become available, it may be worthwhile to investigate the type of data assimilation scheme that would minimize wind speed errors.Furthermore, higher resolution reanalysis need not necessarily be superior for initial and boundary conditions and it is possible that data assimilation schemes used in producing the coarse input atmospheric data may be a deciding factor even for high resolution wind speed predictions.However, the wind speed error distributions differentiate between large hourly negative and positive errors which are averaged to determine the MBE statistic.The error distribution/histogram determines how often the NWP make unacceptable errors.Model results at Crown Point with NCEP as input data produced 86% of modeled hourly values within ±2 m/s of the actual values while those with FNL input data produced 64% (Figure 12).The percentage of hourly modeled wind speeds within ±1 m/s of the actual values were 55% with NCEP input data and 34% for FNL input data.At Piarco, WRF results using either FNL or NCEP as input have approximately 87% of modeled hourly values within ±2 m/s of the actual wind speed and 56% within ±1 m/s (Figure 13).
In summary, the coarser NCEP-DOE reanalysis (1.875° horizontal resolution) provides more accurate estimates of the wind speeds than the FNL (1° resolution) analysis in terms of minimizing mean wind speed errors and producing more symmetric wind speed error distributions for the two measurement stations in Trinidad and Tobago.Sensitivity tests for input and boundary conditions   However, the wind speed error distributions differentiate between large hourly negative and positive errors which are averaged to determine the MBE statistic.The error distribution/histogram determines how often the NWP make unacceptable errors.Model results at Crown Point with NCEP as input data produced 86% of modeled hourly values within ±2 m/s of the actual values while those with FNL input data produced 64% (Figure 12).The percentage of hourly modeled wind speeds within ±1 m/s of the actual values were 55% with NCEP input data and 34% for FNL input data.At Piarco, WRF results using either FNL or NCEP as input have approximately 87% of modeled hourly values within ±2 m/s of the actual wind speed and 56% within ±1 m/s (Figure 13).
In summary, the coarser NCEP-DOE reanalysis (1.875° horizontal resolution) provides more accurate estimates of the wind speeds than the FNL (1° resolution) analysis in terms of minimizing mean wind speed errors and producing more symmetric wind speed error distributions for the two measurement stations in Trinidad and Tobago.Sensitivity tests for input and boundary conditions

Sensitivity to PBL Schemes
WRF was run for seven PBL schemes for the month of April 2003 on a grid of 5 km horizontal resolution which is nested within a 25 km resolution parent grid.The mean error statistics in wind speed and wind direction at each site, Table 2 for Piarco and Table 3 for Crown Point, are based on hourly averaged wind speed data from the only two WMO meteorological stations in the southernmost Caribbean islands of Trinidad and Tobago.The comparisons are performed for the 5 km horizontal resolution and we expect that comparative results among the various PBL schemes to hold at the 1 km horizontal resolution as preliminary tests indicate that there is little improvement in error statistics in increasing horizontal resolution from 5 km to 1 km.The results focus on wind speeds rather than wind directions because the main application of the sensitivity tests of this work is in mapping the wind resources of the Trinidad and Tobago.Tables 2 and 3 have two significant features: firstly, against most wind speed error statistics, the TEMF PBL scheme performed the worst at the two sites, and secondly, with the exception of TEMF, no one PBL scheme seems to perform better than any other PBL scheme.The TEMF scheme produced the largest RMSE, MABE, STDE, second largest MBE, smallest R 2 coefficient and IOA in wind speed.With the exception of TEMF, the RMSE, MABE, and STDE in wind speed for the other PBL schemes lie within a small range: the RMSE in wind speed was 1.28-1.47m/s at Piarco and 1.76-2.11m/s at Crown Point, the MABE was 1.02-1.16m/s at Piarco and 1.39-1.52m/s at Crown Point, and STDE was 1.26-1.34m/s at Piarco and 1.67-1.81m/s at Crown Point.No one PBL scheme produced the smallest bias error in wind speed at both stations.MYNN 2.5 and UW-TK both produced the smallest RMSE at Piarco and ACM2 had the smallest RMSE at Crown Point.The smallest MABE at Piarco was from the MYNN 2.5 PBL scheme while the smallest MABE at Crown Point resulted from the use of YSU-TOPO.YSU and UW-TK produced the highest R 2 coefficient and IOA at Piarco.QNSE gave the highest R 2 at Crown Point and UW-TK, the highest IOA at Crown Point.As no PBL outperforms the other at each station, the mean errors were calculated, which are shown in Table 4.
YSU-TOPO provides the smallest overestimation of approximately 0.16 m/s with MYNN 2.5 giving the smallest underestimation of about −0.16 m/s.The RMSE is smallest for ACM2, YSU-NO-TOPO and YSU-TOPO.These PBL schemes also produced the smallest MABE.With the exception of TEMF, the STDE are similar for all PBL schemes tested with PBLs YSU-TOPO, ACM2, UW-TK, and QNSE having approximately the smallest values.The average spatial error statistics therefore indicate that the YSU-TOPO scheme gives the smallest overall spatial errors.PBL schemes were also assessed in their abilities to model the diurnal variation in wind speed.Most PBL schemes were able to replicate the smaller diurnal amplitude at Crown Point (Figure 14) relative to Piarco (Figure 15) with the exception of TEMF which significantly underestimates the diurnal amplitude at both sites.Although the YSU-TOPO was chosen as the PBL scheme with best overall spatial performance, it should be noted that YSU overestimates the night-time wind speeds from 10 p.m. to 7 a.m. at both sites with small positive biases from 8 a.m. to 5 p.m. (Figures 14 and 15).
In addition, YSU-NO-TOPO and YSU-TOPO produced similar magnitudes in diurnal variation at both sites and along with MYNN 2.5, UW-TK, QNSE and ACM2, these schemes generated daytime maximum wind speeds similar to the observed maximum.For the two sites considered QNSE, TEMF, and UW-TK gave similar nighttime wind speeds and have smaller biases for nighttime winds than other PBL schemes.The diurnal variation at the Piarco site for all PBLs (Figure 15) show large positive biases for wind speeds less than 2 m/s which occur at night.The PBL schemes do not capture the mean nighttime wind speeds.This has been observed in several other studies [16,17] and may be linked to the stable conditions and low turbulence that persist during the night but with the PBL schemes producing too much vertical mixing at night [16].The PBL ensemble, however, captures the hourly wind speeds between 8 a.m. and 7 p.m. at Crown Point and 9 a.m. and 5 p.m. at Piarco.Various PBL schemes perform better for each hour.For example, at Piarco, UW-TK best represented the mean hourly wind speeds from 10 p.m.to a.m., MYNN 2.5 from 8 a.m. to 1 p.m., YSU-TOPO from 2 p.m. to 4 p.m., MYJ from 5 p.m.to 6 p.m., and TEMF from 7 p.m. to 9 p.m. Table 4. Mean spatial error statistics in wind speed at the 10 m height above ground for each PBL scheme.These spatial errors were derived from averaging the corresponding station errors in Tables 3  and 4. QNSE having approximately the smallest values.The average spatial error statistics therefore indicate that the YSU-TOPO scheme gives the smallest overall spatial errors.PBL schemes were also assessed in their abilities to model the diurnal variation in wind speed.Most PBL schemes were able to replicate the smaller diurnal amplitude at Crown Point (Figure 14) relative to Piarco (Figure 15) with the exception of TEMF which significantly underestimates the diurnal amplitude at both sites.

Statistic
Although the YSU-TOPO was chosen as the PBL scheme with best overall spatial performance, it should be noted that YSU overestimates the night-time wind speeds from 10 p.m.-7 a.m. at both sites with small positive biases from 8 a.m. to 5 p.m. (Figures 14 and 15).In addition, YSU-NO-TOPO and YSU-TOPO produced similar magnitudes in diurnal variation at both sites and along with MYNN 2.5, UW-TK, QNSE and ACM2, these schemes generated daytime maximum wind speeds similar to the observed maximum.For the two sites considered QNSE, TEMF, and UW-TK gave similar nighttime wind speeds and have smaller biases for nighttime winds than other PBL schemes.The diurnal variation at the Piarco site for all PBLs (Figure 15) show large positive biases for wind speeds less than 2 m/s which occur at night.The PBL schemes do not capture the mean nighttime wind speeds.This has been observed in several other studies [16,17] and may be linked to the stable conditions and low turbulence that persist during the night but with the PBL schemes producing too much vertical mixing at night    An accurate prediction of the wind speeds is crucial in determining the wind power density (WPD) of each site.The WPD is taken as a measure of the wind potential or resource of a site.Sites with higher wind potential are of greater interest to developers.In this study, the WPD was determined using the standard air density to compare the impact of the PBL-SL predictions of wind speed on WPD.Although the YSU PBL scheme was found to have a positive bias in the wind speed, it provides a conservative estimate and reasonable magnitude of the WPD at the two measurement sites.The YSU-TOPO PBL scheme underestimates the WPD by 2.5 W m −2 at Piarco, even though it predicts a small positive bias in mean monthly wind speed.Although the MYNN 2.5 provides a more conservative estimate for mean wind speed over YSU, it underestimates the mean WPD by 22.1 W m −2 .At Crown Point, MYNN 2.5 underestimates the WPD by 32.9 W m −2 compared with YSU's underestimate of 18.0 W m −2 .In contrast, ACM2 overestimates the mean WPD by 4.0 W m −2 at Piarco and 2.9 W m −2 at Crown Point.
In addition to accurately predicting the diurnal evolution of wind speeds, the selected PBL model has to be able to simulate wind speed histograms as wind speed distributions directly impact on the estimate of the WPD.The wind speed distributions are often approximated by probability frequency distributions to estimate the energy output of wind turbines.The two sites, Crown Point and Piarco, are characterized by bimodal wind speed distributions as denoted by the red lines on Figures 16 and 17 An accurate prediction of the wind speeds is crucial in determining the wind power density (WPD) of each site.The WPD is taken as a measure of the wind potential or resource of a site.Sites with higher wind potential are of greater interest to developers.In this study, the WPD was determined using the standard air density to compare the impact of the PBL-SL predictions of wind speed on WPD.Although the YSU PBL scheme was found to have a positive bias in the wind speed, it provides a conservative estimate and reasonable magnitude of the WPD at the two measurement sites.The YSU-TOPO PBL scheme underestimates the WPD by 2.5 W m −2 at Piarco, even though it predicts a small positive bias in mean monthly wind speed.Although the MYNN 2.5 provides a more conservative estimate for mean wind speed over YSU, it underestimates the mean WPD by 22.1 W m −2 .At Crown Point, MYNN 2.5 underestimates the WPD by 32.9 W m −2 compared with YSU's underestimate of 18.0 W m −2 .In contrast, ACM2 overestimates the mean WPD by 4.0 W m −2 at Piarco and 2.9 W m −2 at Crown Point.
In addition to accurately predicting the diurnal evolution of wind speeds, the selected PBL model has to be able to simulate wind speed histograms as wind speed distributions directly impact on the estimate of the WPD.The wind speed distributions are often approximated by probability frequency distributions to estimate the energy output of wind turbines.The two sites, Crown Point and Piarco, are characterized by bimodal wind speed distributions as denoted by the red lines on Figures 16 and 17.The TEMF and QNSE PBL schemes do not capture the bimodal distribution at Crown Point and Piarco.The MYNN 2.5 also shows such undesired behavior at Crown Point.At Crown Point, YSU-NO-TOPO, YSU-TOPO, and ACM2 underestimate the frequency of the low wind speed intervals 0-1 m/s and 1-2 m/s and overestimate the 5-6 m/s interval at the Crown Point site (Figure 16).The MYJ better represents the 2-3 m/s and 5-6 m/s intervals.MYJ, YSU, YSU-TOPO, ACM2 captures the wind speed intervals in which the peaks occur.UW-TK overestimates the frequency of wind speed bins less than 4 m/s at Crown Point.At Piarco, ACM2 and UW-TK overestimates the 2-3 m/s interval and underestimates the 4-5 m/s interval frequency.Of the remaining PBL schemes, the MYJ PBL scheme captures the intervals that give the observed frequency distribution its bimodal nature, but underestimates the occurrence of the 0-1 m/s and 1-2 m/s intervals (Figure 17).The YSU-NO-TOPO, YSU-TOPO, and the MYJ accurately estimate the frequency of the 6-7 m/s wind speed interval, which constitutes the second peak in the wind speed frequency distribution.MYJ overestimates the frequency of the 7-8 m/s while the YSU schemes estimate its frequency more accurately.The YSU schemes tend to overestimate the frequency of the 2-3 m/s, 3-4 m/s, and 5-6 m/s interval frequencies and underestimate the frequency of the 4-5 m/s interval.
Energies 2017, 10, 931 18 of 24 distribution its bimodal nature, but underestimates the occurrence of the 0-1 m/s and 1-2 m/s intervals (Figure 17).The YSU-NO-TOPO, YSU-TOPO, and the MYJ accurately estimate the frequency of the 6-7 m/s wind speed interval, which constitutes the second peak in the wind speed frequency distribution.MYJ overestimates the frequency of the 7-8 m/s while the YSU schemes estimate its frequency more accurately.The YSU schemes tend to overestimate the frequency of the 2-3 m/s, 3-4 m/s, and 5-6 m/s interval frequencies and underestimate the frequency of the 4-5 m/s interval.were able to capture the bimodal feature of the wind speed distribution.The bimodal feature is due to the high occurrence of low wind speeds at the 10 m level.Although this qualitative assessment was performed for 10 m level wind speed distributions, it is here hypothesized that schemes capable of capturing the bimodal wind speed distribution at the surface would most likely be able to capture the typical unimodal wind speeds at higher levels above ground.

Conclusions
In this study, various model configurations of the WRF model that minimize the errors in the wind resources of the southernmost Caribbean islands of Trinidad and Tobago were investigated.The general conclusion is that the model has the ability to capture some of the essential features of diurnal variability in wind speed, wind speed histograms and wind direction distributions.Specific aspects of the model that were optimized were the spin-up period, the input atmospheric conditions, and the planetary boundary layer (PBL) scheme.The WRF configuration of a 12 h spin-up and the YSU-TOPO PBL parameterization scheme minimize the spatial average errors when used in conjunction with the RRTM scheme for longwave radiation, the Dudhia scheme for shortwave radiation, the thermal diffusion land surface model, the WSM 3 scheme for microphysics on the inner nests and the Betts-Miller-Janjic scheme for cumulus parameterization on the outermost parent domain.In addition, the NCEP reanalysis was found to be more suitable as input initial and boundary conditions over the FNL analyses.This configuration will be used in subsequent studies for developing numerical wind maps for the Caribbean islands of Trinidad and Tobago.It is noted here that the observational data sets used to optimize the WRF model were limited to the 10 m level and did not include wind data for greater heights more relevant to commercial wind turbines.However, as the model configuration is able to capture a complex feature of the 10 m wind speed histogram, that is, the bimodal distribution characteristic, the suggested model configuration may be able to capture the typical unimodal wind speed distributions at greater heights.Future tall-tower experiments will be required to supplement the 10 m wind data currently available for the southernmost Caribbean islands of Trinidad and Tobago to verify whether the model configuration of the WRF model optimized for the 10 m level wind resources would be able to optimize the wind resources at greater heights.

24 Figure 1 .
Figure 1.Numerical domain under consideration.The outermost boundaries denote the parent grid, d02 and d03 are the nested domains, d03 encompasses the area of interest, Trinidad and Tobago.

Figure 2 .
Figure 2. Computational grid in the vertical.(a) Profile of all vertical grid levels through a cross-section across 10.72° N along which the highest topographic elevation (712 m) is located; (b) The lowest nine vertical grid levels as in (a) with the first level indicated in blue and the elevation profile in red; (c) The average height of the numerical grid's vertical levels.

Figure 1 . 24 Figure 1 .
Figure 1.Numerical domain under consideration.The outermost boundaries denote the parent grid, d02 and d03 are the nested domains, d03 encompasses the area of interest, Trinidad and Tobago.

Figure 2 .
Figure 2. Computational grid in the vertical.(a) Profile of all vertical grid levels through a cross-section across 10.72° N along which the highest topographic elevation (712 m) is located; (b) The lowest nine vertical grid levels as in (a) with the first level indicated in blue and the elevation profile in red; (c) The average height of the numerical grid's vertical levels.

Figure 2 .
Figure 2. Computational grid in the vertical.(a) Profile of all vertical grid levels through a cross-section across 10.72 • N along which the highest topographic elevation (712 m) is located; (b) The lowest nine vertical grid levels as in (a) with the first level indicated in blue and the elevation profile in red; (c) The average height of the numerical grid's vertical levels.

Figure 3 .
Figure 3. Kinetic energy spectra at several spin-up times when FNL initial and boundary conditions (ICs/BCs) were used.

Figure 3 .
Figure 3. Kinetic energy spectra at several spin-up times when FNL initial and boundary conditions (ICs/BCs) were used.

Figure 4 .
Figure 4. Kinetic energy spectra at several spin-up times when NCEP ICs/BCs were used.

Figure 5 .
Figure 5. 10 m level wind speed time series at Crown Point for 1-14 April 2003 from in-situ data (red line) compared with WRF model simulations initialized with FNL data (blue line) and those with NCEP-DOE reanalysis data (black line).

Figure 4 .
Figure 4. Kinetic energy spectra at several spin-up times when NCEP ICs/BCs were used.

3. 2 .
FNL vs NCEP As Initial and Boundary Conditions 3.2.1.General Observations for Hourly Simulations of Wind Speed The modeled hourly time series for the first fourteen days of April 2003 are shown in Figures 5 and 6 for Crown Point and Piarco, respectively.The WRF model's configuration with the YSU PBL and with both NCEP and FNL initial conditions (ICs) is capable of reproducing the diurnal behavior of the wind speed at both stations.However, the configurations seem to underestimate the diurnal maxima at Crown Point.The model overestimates the nighttime low wind speeds at both stations.The model underestimates wind speeds that exceed 8 m/s at Piarco.In general, the model configuration and the input atmospheric data capture the timing of the daily maxima and minima moreso at Piarco than Crown Point.There are however, some timing problems at Crown Point on 4 April (~72 h) and 12 April (~280 h).Energies 2017, 10, 931 8 of 24

Figure 4 .
Figure 4. Kinetic energy spectra at several spin-up times when NCEP ICs/BCs were used.

Figure 5 .
Figure 5. 10 m level wind speed time series at Crown Point for 1-14 April 2003 from in-situ data (red line) compared with WRF model simulations initialized with FNL data (blue line) and those with NCEP-DOE reanalysis data (black line).

Figure 5 .
Figure 5. 10 m level wind speed time series at Crown Point for 1-14 April 2003 from in-situ data (red line) compared with WRF model simulations initialized with FNL data (blue line) and those with NCEP-DOE reanalysis data (black line).

Figure 6 .
Figure 6. 10 m level wind speed time series at Piarco for 1-14 April 2003 from in-situ data (red line) compared with WRF model simulations initialized with FNL data (blue line) and those with NCEP-DOE reanalysis data (black line).

Figure 6 .
Figure 6. 10 m level wind speed time series at Piarco for 1-14 April 2003 from in-situ data (red line) compared with WRF model simulations initialized with FNL data (blue line) and those with NCEP-DOE reanalysis data (black line).

Energies 2017, 10 , 931 10 of 24 Figure 7 .
Figure 7. 10 m level wind speed error variations by hour for April 2003 for Crown Point.

Figure 7 .
Figure 7. 10 m level wind speed error variations by hour for April 2003 for Crown Point.

Figure 7 .
10 m level wind speed error variations by hour for April 2003 for Crown Point.

Figure 8 .
Figure 8. 10 m level wind speed error variations by hour for April 2003 for Piarco.

Figure 8 .
Figure 8. 10 m level wind speed error variations by hour for April 2003 for Piarco.

Figure 9 .
Figure 9. 10 m level wind direction distribution at (a) Crown Point and (b) Piarco from in-situ data (red bars) compared with WRF model simulations initialized with FNL data (blue bars) and those with NCEP reanalysis data (black bars).

Figure 9 .
Figure 9. 10 m level wind direction distribution at (a) Crown Point and (b) Piarco from in-situ data (red bars) compared with WRF model simulations initialized with FNL data (blue bars) and those with NCEP reanalysis data (black bars).

Figure 10 .
Figure 10.Histograms of 10 m level wind speed at Crown Point as simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 11 .
Figure 11.Histograms of 10 m level wind speed at Piarco simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 10 .
Figure 10.Histograms of 10 m level wind speed at Crown Point as simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 10 .
Figure 10.Histograms of 10 m level wind speed at Crown Point as simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 11 .
Figure 11.Histograms of 10 m level wind speed at Piarco simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 11 .
Figure 11.Histograms of 10 m level wind speed at Piarco simulated by WRF using FNL as input atmospheric conditions (in blue), with NCEP as input (in black) as compared against in-situ observations (in red).

Figure 12 .
Figure 12.Wind speed error distributions at Crown Point for the 10 m height above ground.

Figure 13 .
Figure 13.Wind speed error distributions at Piarco for the 10 m height above ground.

Figure 12 .
Figure 12.Wind speed error distributions at Crown Point for the 10 m height above ground.

Figure 12 .
Figure 12.Wind speed error distributions at Crown Point for the 10 m height above ground.

Figure 13 .
Figure 13.Wind speed error distributions at Piarco for the 10 m height above ground.

Figure 13 .
Figure 13.Wind speed error distributions at Piarco for the 10 m height above ground.

Figure 14 .
Figure 14.Diurnal variation in 10 m level wind speed as simulated by 7 PBL schemes and compared with in-situ observations at Crown Point.

Figure 14 .
Figure 14.Diurnal variation in 10 m level wind speed as simulated by 7 PBL schemes and compared with in-situ observations at Crown Point.

Figure 15 .
Figure 15.Diurnal variation in 10 m level wind speed as simulated by 7 PBL schemes and compared with in-situ observations at Piarco.

Figure 15 .
Figure 15.Diurnal variation in 10 m level wind speed as simulated by 7 PBL schemes and compared with in-situ observations at Piarco.

Figure 16 .
Figure 16.Histograms of simulated 10 m level wind speeds at Crown Point for April 2003.The red line denotes the observed wind speed histogram for April 2003.

Figure 16 .
Figure 16.Histograms of simulated 10 m level wind speeds at Crown Point for April 2003.The red line denotes the observed wind speed histogram for April 2003.

Table 2 .
Error statistics for wind speed and wind direction at the 10 m height above ground, and an estimate of the monthly wind power density for April 2003 at Piarco.Error statistics were computed for 720 h of simulated winds.

Table 3 .
Error statistics for wind speed and wind direction at the 10 m height above ground, and an estimate of the monthly wind power density for April 2003 at Crown Point.Error statistics were computed for 720 h of simulated winds.