Fig. (1) shows a brief overview of the methodology. The overall idea is to conduct discrete LES that represents day-to-day variations in atmospheric stability using realistic boundary conditions from field data. The LES tool requires discrete values for surface heat flux on the lower boundary, initial inversion height (specified using a 3-layer temperature profile as explained in section 2.2), and geostrophic wind on the upper boundary. The methodology consists of following steps, which will be explained in detail in this section: (1) Obtaining the instantaneous surface heat flux from the field data, (2) gathering the surface heat flux data occurring between 12-13 LT (for each month), (3) fitting the gathered data to multiple Probability Density Functions (PDFs), (4) selecting the best fit PDF using, and (5) determining the discrete boundary conditions required for LES using the surface heat flux PDF. The inversion height for each surface heat flux is determined from the field data correlation between surface heat flux and boundary height as shown in Fig. (5). The geostrophic wind for each surface heat flux is determined from the departure function defined by Clarke and Hess [26Clarke RH, Hess GD. Geostrophic departure and the functions A and B of Rossby-number similarity theory. Boundary-Layer Meteorol 1974; 7: 267-87.
[http://dx.doi.org/10.1007/BF00240832] ].

The PDF for the data occurring between 12LT-13LT represents the day-to-day variations at 12LT for an entire month. The study uses the PDF to determine discrete values for the surface heat flux data occurring between 12LT-13LT. The discrete surface heat flux values become inputs into LES in addition to the corresponding initial inversion height, and geostrophic wind. The LES profiles are compared to the field data wind speed and TKE profiles to validate the methodology. The Correlations of surface heat flux and: TKE, friction velocity, wind shear, and wind speed for the LES simulations and measured data are compared to understand the implications of day-to-day variations on wind turbines.

The data used in this study is from the 135-m tall M5 meteorological tower installed at the NREL NWTC near Boulder Colorado [27 NWTC Information Portal (NWTC 135-m Meteorological Towers Data Repository). https:// nwtc.nrel.gov/ 135mData Last modified 01-April-2015 ; Accessed 21-October-2018.]. The data consists of high-frequency measurements of wind speed, wind direction and air temperature at elevations of 15, 30, 50, 76, 100, and 131 above the ground January, May, July and October of the year 2016. These months correspond to Winter, Spring, Summer, and Fall for that year. For our analysis, the measured data were used to obtain the derived atmospheric variables, namely, Surface heat flux, atmospheric boundary layer height, wind shear, and Turbulent Kinetic Energy (TKE) (see Appendix A for details [24Clifton A, Schreck S, Scott G, Kelley N, Lundquist JK. Turbine inflow characterization at the national wind technology center. J Sol Energy Eng 2013; 135: 031017-031017-11.
[http://dx.doi.org/10.1115/1.4024068] ].

Fig. (2a) shows the monthly-averaged surface heat flux for each hour of the day for January, May, July, and October. The vertical bars denote one standard deviation (+1σ to -1σ) from the mean. Fig. (2b) shows the standard deviation (σ) of the surface heat flux at each hour. In general, a positive value of surface heat flux refers to convective stability state of the atmosphere, whereas negative values denote a stable atmosphere, and a zero value refers to a neutral atmosphere [3Stull RB. An introduction to boundary layer meteorology 2012; Vol. 13]. The strength of the stability is related to the magnitude of the value. Significant differences in both the values of the mean surface flux and the standard deviation from the mean is observed between the 4 selected months. The mean surface heat flux during the early morning hours (00 to 05 LT) is negative for all the 4 months with values ranging from -0.0295 kms^-1 for January to -0.158 kms^-1 for May. Positive mean surface heat flux occurs between times of 9-15 LT, 6-16 LT, 6-17 LT, and 7-15 LT for January, May, July, and October respectively. The maximum mean surface heat flux occurs near 12 LT for all the months, but the magnitude ranges from 0.088, 0.234, 0.231, and 0.177 kms^-1 for January, May, July, and October respectively. During the later hours, the mean surface heat flux transitions back to negative values.

Fig. (1) Methodology to incorporate boundary and initial conditions into LES.

Fig. (2) Mean and variability of surface heat flux at different hours for different months using NREL data. (a) shows the mean and one standard deviation error bars (b) shows the standard deviation.

Fig. (2b) shows the month of October to have the highest standard deviations from the mean at 10 LT with a standard deviation of 0.126 kms^-1 whereas, for July, the highest standard deviation is 0.125 kms^-1 at 12 LT. The highest standard deviation for January and May is 0.070 kms^-1 and 0.124 kms^-1 which occur at 10LT and 12LT respectively. The July 12 LT standard deviation is nearly 58% of the mean surface heat flux at that hour, which indicates the day-to-day variations of the surface flux are significant. It is interesting to observe that though the general trend of mean surface heat flux is similar for all the months; the standard deviation from the mean does not follow the same trend. The difference is due to the fact the maximum mean surface heat flux for the summer month of July is over two times the value than January. The difference in maximum mean surface heat flux occurs even though the minimum values of the mean surface heat flux are comparable between all the 4 months. July exhibits the highest peak surface flux, but the variability is highest for May, July, and May.

The magnitude and variation of surface heat flux changes between different seasons, especially during the convective portion of the day from 8 to 16 LT. July shows the largest mean surface heat flux of 0.24 kms^-1 with its highest standard deviations of 0.125 kms^-1 at 12 LT while January shows the lowest mean surface heat flux of 0.088 kms^-1 and the lowest standard deviation of 0.06 kms^-1 at 12 LT. The results show clear seasonal (differences between the four months) and diurnal variations in surface heat flux exist between these four months. In addition, the surface heat flux standard deviation from the mean shows a day-to-day variation as well.

The day-to-day variations over multiple diurnal cycles are studied. For this purpose, a scatter plot of 10-min averaged surface heat flux for the entire month at a given time is analyzed. Fig. (3a) shows the 10-minute averaged surface heat flux for averaged data that falls between 12 LT and 13 LT time for all the days of July (referred to as July 12LT). The x-axis of is the day of the month for July 2016. The figure shows significant day-to-day variations in surface heat flux with the highest value of 0.52 kms^-1 and the lowest value of 0.05 kms^-1. The day-to-day variations observed in the average surface heat flux for the July 12 LT can be well represented using a Probability Density Function (PDF) using mathematical probability theory [28Hogg RV, Craig AT. Introduction to mathematical statistics(5th edition) 1995.]. A normal PDF fitted to the data from July 12 LT is shown in Fig. (3b). Fitted in this manner, 68% of the data is within the ±1σ, 95% of the data is within ±2σ, and 99% fall within ±3σ.

The PDF represents the variability at 12LT and will be used to select the inputs for LES. The actual PDF is determined in section 3.1. The focus of the current study is on day-to-day variations, particularly for the July 12 LT, because of the large mean surface heat flux and standard deviation from the mean that occurs at that time. In addition, the January 12 LT will also be studied to investigate the implications of smaller mean and standard deviation from the mean of the surface heat flux.

Fig. (3) (a) 10 min averaged surface heat flux for the July 12 LT local time for each day (b) probability density function for the July 12 LT local time surface heat flux.

Next, we investigate the correlations of the surface heat flux with TKE and ABL height. The correlation with TKE helps to understand how the variations in surface heat flux will impact the output of the LES. Fig. (4) shows the correlation of 10-min averaged surface heat flux and TKE for different months of the NREL data at an elevation of 76 m above the ground. There is a positive correlation between surface heat flux and TKE (R² of 0.45). The low R² indicates some variability, but the correlation has a p-value <0.01 indicating the correlation is significant even with the variability. The methodology to incorporate surface heat flux variations into LES will need to capture the relationships shown to be used in LES applications. The surface heat flux correlations will be used to validate the LES in section 3.3.

Fig. (4) Correlation of heat flux and TKE at 76 m during convective hours for different months.

Finally, we investigate the correlation between surface heat flux and the approximated ABL height (shown in Eq. 16). The relationship is shown by plotting the two on a log-log scaled plot as seen in Fig. (5). The data show a positive correlation. The correlation function between surface heat flux and ABL height is a power law, described as , rather than linear as found with TKE. This relationship had slightly less variability than TKE with an R² value of 0.55 compared to 0.47. Again, the p-value is <0.01, which indicates significance to the correlation. The correlation function is used to obtain the initial inversion height for LES simulations based on the boundary condition given for surface heat flux [29Sugiyama G, Nasstrom JS. Methods for determining the height of the atmospheric boundary layer. Lawrence Livermore Natl. Lab Intern Rep UCRL-ID-133200 1999.].

Fig. (5) Correlation of surface heat flux and ABL height during convective hours for different months.

In summary, the NWTC Metrological tower data for January, May, October, and December 2016 shows the presence of day-to-day variability of surface heat flux over multiple diurnal cycles. Results have been shown for July 12 LT. However, all four months show similar trends (section 3.1). These results are comparable to Grimmond and Oke [11Grimmond CSB, Oke TR. Comparison of heat fluxes from summertime observations in the suburbs of four North American cities. J Appl Meteorol 1995; 34: 873-89.
[http://dx.doi.org/10.1175/1520-0450(1995)034<0873:COHFFS>2.0.CO;2] ] who showed the largest day-to-day variations of surface heat flux occur near midday. The 10-minute averaged surface heat flux showed correlations with both TKE and ABL height, which is also described by Obukhov theory [30Foken T. 50 years of the Monin–Obukhov similarity theory. Boundary-Layer Meteorol 2006; 119: 431-47.
[http://dx.doi.org/10.1007/s10546-006-9048-6] ]. A linear correlation exists between the averaged surface heat flux and TKE, whereas the correlation is a power law with the ABL height. The correlations have implications to the modeling of LES. The variations of surface heat flux affect the TKE output of the simulation, which is important in may LES applications. The study investigates multiple PDFs to best define the variations in surface heat flux. The study uses multiple discrete LES simulations determined from a mathematical PDF (section 3.1) that represent the monthly day-to-day variations, rather than a one mean value. The LES methodology is described in subsection 2.2.

The Large Eddy Simulation (LES) has been performed using the Simulator for Wind Farm Applications (SOWFA) solver developed by the National Renewable Energy Laboratory (NREL). The LES uses ideal assumptions with a homogenous pressure gradient, flat terrain, and steady conditions (Churchfield et al. [13Churchfield M J, Moriarty P J, Vijayakumar G, Brasseur J G. Wind Energy-Related Atmospheric Boundary-Layer Large-Eddy Simulation Using OpenFOAM. Natl. Renew Energy Lab Gold CO Rep. No NRELCP-500-48905 2010., 31Churchfield MJ, Lee S, Moriarty PJ, et al. A large-eddy simulation of wind-plant aerodynamics. AIAA Pap 2012; 537: 2012.]). For the sake of completeness, the details of the solver are repeated here. The governing equations are the filtered Navier-Stokes (momentum and potential temperature transport equations) with Boussinesq approximations as follows:

(1)

(2)

(3)

where, a tilde denotes the spatially resolved filtered component at scale, denotes the modified pressure. The spatial gradient of the mean pressure term acts to drive the flow convection, (u₁, u₂, u₃) = (u, v, w) are the components of the velocity field (u is the axial component; v is the lateral or y-component, and w is the vertical component). θ is the resolved potential temperature; θ_O is the reference temperature. The Subgrid Stresses (SGS) are included in the τ term along with the fluid stress tensor. The SGS flux and SFS temperature flux are parameterized, respectively as follows:

where D is the deviatoric SGS stress tensor; v^SGS is the SGS eddy viscosity; and Pr^SGS is the SGS Prandtl number. The SGS flux is computed using the Moeng’s model [14Ayotte KW, Sullivan PP, Andrén A, et al. An evaluation of neutral and convective planetary boundary-layer parameterizations relative to large eddy simulations. Boundary-Layer Meteorol 1996; 79: 131-75.
[http://dx.doi.org/10.1007/BF00120078] ], where the turbulent viscosity and turbulent kinetic energy are modeled as follows:

using the mesh cells in x, y and z directions (the resolved strain-rate tensor is) and Pr_SGS is constant (ref dependent on stratification level).

The finite volume solver uses the Pressure Implicit with Splitting of Operations (PISO) algorithm to solve the momentum equation [32Issa R. Solution of the implicit discretized fluid flow equations by operator splitting (PDF Download Available) Available online: https:// www.researchgate.net/ publication/ 247576431_ Solution_of _the_ Implicit_ Discretized_ Fluid _Flow _Equations _by _Operator_ Splitting (Accessed on Dec 29, 2016).]. The pressure and velocity are solved implicitly in this method using Crank-Nicholson time advancement for the predictor and two corrector steps [13Churchfield M J, Moriarty P J, Vijayakumar G, Brasseur J G. Wind Energy-Related Atmospheric Boundary-Layer Large-Eddy Simulation Using OpenFOAM. Natl. Renew Energy Lab Gold CO Rep. No NRELCP-500-48905 2010.]. The code solves the buoyant term, the SGS viscosity, and the Coriolis terms explicitly. The pressure gradient, , drives the flow and is constant in space. The pressure gradient is adjusted to maintain the geostrophic wind at the top of the domain. The Coriolis parameter, f^c, is calculated based on the latitude, which is 40° where the NREL site is located [17Kumar V, Svensson G, Holtslag AM, Meneveau C, Parlange MB. Impact of surface flux formulations and geostrophic forcing on large-eddy simulations of diurnal atmospheric boundary layer flow. J Appl Meteorol Climatol 2010; 49: 1496-516.
[http://dx.doi.org/10.1175/2010JAMC2145.1] ]. The deviatoric SGS stress tensor, , is equal to where v^SGS is the SGS viscosity and is the resolved strain rate tensor. The SGS viscosity is modeled using the one-equation SGS model given by Moeng [4Moeng C-H. A large-eddy-simulation model for the study of planetary boundary-layer turbulence. J Atmos Sci 1984; 41: 2052-62.
[http://dx.doi.org/10.1175/1520-0469(1984)041<2052:ALESMF>2.0.CO;2] ].

Each simulation is initialized with an initial temperature and velocity profile shown in Eq. 4. where C is a constant and is 280 for the January simulations and 300 for the July simulations. The initial temperature profile is a three-layer structure [33Sullivan PP, Patton EG. The effect of mesh resolution on convective boundary layer statistics and structures generated by large-eddy simulation. J Atmos Sci 2011; 68: 2395-415.
[http://dx.doi.org/10.1175/JAS-D-10-05010.1] ]. The first is a constant temperature gradient to the inversion height. At the initial inversion height, the temperature is changed by 8° K over 150 m. Finally, a constant temperature gradient of 0.003 Km^-1 is given from the top of the inversion height to the top of the domain [13Churchfield M J, Moriarty P J, Vijayakumar G, Brasseur J G. Wind Energy-Related Atmospheric Boundary-Layer Large-Eddy Simulation Using OpenFOAM. Natl. Renew Energy Lab Gold CO Rep. No NRELCP-500-48905 2010.]. The initial velocity is prescribed as a constant geostrophic wind velocity throughout the domain. Random temperature perturbations and a sinusoidal velocity perturbation are added near the surface to decrease the amount of time to generate realistic turbulence. As the simulation progresses, the profiles develop to those of an ABL.

(4)

All lateral boundary conditions are specified as periodic. The upper-velocity boundary condition is specified as zero velocity normal to the boundary and zero gradient parallel to the boundary. A Neumann boundary is prescribed as the upper-pressure boundary condition that is determined from the normal momentum equation. The upper-temperature boundary condition is specified as a constant temperature gradient of 0.003 to match the initialized temperature profile. The upper SGS stress boundary condition is given as a zero gradient.

The surface roughness z_O [8Paulson CA. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J Appl Meteorol 1970; 9: 857-61.
[http://dx.doi.org/10.1175/1520-0450(1970)009<0857:TMROWS>2.0.CO;2] ], which represents different land and surface types, sets the lower velocity boundary condition. z_O is given as 0.15 m The lower pressure boundary condition is given as the gradient determined by the normal momentum equation. The lower temperature boundary condition is given as a constant temperature flux (surface heat flux). The lower SGS viscosity boundary condition is determined using the Schumann Model [34Schumann U. Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. J Comput Phys 1975; 18: 376-404.
[http://dx.doi.org/10.1016/0021-9991(75)90093-5] ]. The SGS viscosity values are determined to satisfy the Schumann Model surface stress as described by Churchfield et al. [13Churchfield M J, Moriarty P J, Vijayakumar G, Brasseur J G. Wind Energy-Related Atmospheric Boundary-Layer Large-Eddy Simulation Using OpenFOAM. Natl. Renew Energy Lab Gold CO Rep. No NRELCP-500-48905 2010.].

The domain of the simulation is 4000 m x 4000 m x 2000 m. The predominant wind is set at an angle through the domain such that the large-scale structures can properly form [13Churchfield M J, Moriarty P J, Vijayakumar G, Brasseur J G. Wind Energy-Related Atmospheric Boundary-Layer Large-Eddy Simulation Using OpenFOAM. Natl. Renew Energy Lab Gold CO Rep. No NRELCP-500-48905 2010.] and the turbulent statistics did not change with increased domain size. Simulations were performed with grid points of Nx, Ny, and Nz as 160, 160 and 128 in the horizontal (x), spanwise (y), and vertical (z) - directions respectively. The grid points in the horizontal direction are equally spaced while the spacing in the vertical direction increases such that the last cell is four times larger than the first cell where Δz = .146 Nz + 5.98. The simulations were run for 9,000s and then the statistics were taken for another 2,000s (the eddy turnover times ranged from 6-12 min). Table 1 gives a list of LES cases performed for the validation. These were determined using the PDFs for surface heat flux as described in Section 3.5.

Table 1
The LES cases used for validation which were selected by using the July and January PDFs acquired from field data. The table shows the boundary and initial conditions for each case.

The Probability Density Functions (PDFs) are generated by grouping together surface heat flux values over a given hour, applying multiple types of PDFs to the set of data, and selecting the PDF with the best fit as determined by an R² value [35Morgan EC, Lackner M, Vogel RM, Baise LG. Probability distributions for offshore wind speeds. Energy Convers Manage 2011; 52: 15-26.
[http://dx.doi.org/10.1016/j.enconman.2010.06.015] ]. Normal and logistic distributions were tested. The equations for the PDFs are:

(5)

(6)

where µ is the mean, σ is the standard deviation, and s is the shape parameter.

One method to validate the fits of the PDFs is by using the cumulative density function. The cumulative density function represents the probability a random variable X will be less than or equal to x. The field data is compared to the PDF, and an R² can be determined. Fig. (6) shows the cumulative density function fits for the PDFs of surface heat flux. Fig. (6) shows the cumulative density function fit for the (a) Normal PDF and (b) Logistic PDF using data from July at hour 12 LT. Fig. (6) also shows the cumulative density function fit for the (c) Normal PDF and (d) Logistic PDF data from July at hour 21 LT. The figure shows the 12 LT and 21 LT as representative examples for both a convective and a stable instance in the ABL. Most of the deviation between the fit and the observation occurs near the extremes of the cumulative density function. Appendix B gives a complete list of R² values for each hour of each month. Tables 4-7 record the probability functions and their parameters at each hour for the January, May, July, and October data sets. The table shows that the normal distribution had a better fit during the convective portion of the day between 8 and 16 LT. The convective part of the day also had higher R² values ranging from 0.93 to 0.99. The logistic distribution had a better fit during night time conditions between 18 and 6 LT. During this time the R² values were much smaller ranging from 0.61 to 0.89 (which represents a poor fit). The convective ABL data showed a strong fit with the normal distribution. Therefore, the normal PDF was used to select discrete inputs for LES.

Fig. (6) Cumulative probability with fits (a) NREL July hour 14 LT Normal Distribution (b) NREL July hour 14 LT Logistic Distribution (c) NREL July hour 21 LT Normal Distribution (d) NREL July hour 21 LT Logistic Distribution.

The study simulated LES cases for July 12 LT and January 12 LT using the discrete surface flux value obtained from the normal PDF. The simulations consisted of the mean, -1 standard deviation and +1 standard deviation surface heat flux for each July 12 LT and January 12 LT normal PDF. For July 12 LT the mean surface heat flux mean was 0.23 kms^-1 with a standard deviation of 0.125 kms^-1. The mean surface flux in January 12 LT is 0.088 kms^-1 with a standard deviation was 0.06 kms^-1. The July 12LT was selected because it had the largest mean surface heat flux and near the largest standard deviation. The January 12 LT was selected for comparison because the mean and standard deviation were much smaller. The correlation between boundary layer height and surface heat flux shown in Fig. (5) was then used to calculate the boundary layer height for each of the six surface heat fluxes. The boundary layer heights set the initial inversion height. Finally, the geostrophic wind was estimated using u^* in the geostrophic departure described by Clarke and Hess [26Clarke RH, Hess GD. Geostrophic departure and the functions A and B of Rossby-number similarity theory. Boundary-Layer Meteorol 1974; 7: 267-87.
[http://dx.doi.org/10.1007/BF00240832] ] where the A and B coefficients are experimental values used to correct the departure function. The assumption was made to use the same values found experimentally at noon by Clarke and Hess.

The LES boundary and initial conditions for the six cases are in Table 1. The table shows the case number, the selected month and surface heat flux from the PDF, the surface heat flux used as a boundary condition, the inversion height used as an initial condition, and the geostrophic wind used to drive the flow. The table also shows the corresponding Monin-Obukhov Length from each simulation. Neutral conditions generally occur near -300 m [36Iungo GV, Porté-Agel F. Volumetric lidar scanning of wind turbine wakes under convective and neutral atmospheric stability regimes. J Atmos Ocean Technol 2014; 31: 2035-48.
[http://dx.doi.org/10.1175/JTECH-D-13-00252.1] ]. All values for the Monin-Obukhov Length fall between -20 m and -43 m and are therefore portions of the convective ABL [5Johansson C, Smedman A-S, Högström U, Brasseur JG, Khanna S. Critical test of the validity of Monin–Obukhov similarity during convective conditions. J Atmos Sci 2001; 58: 1549-66.
[http://dx.doi.org/10.1175/1520-0469(2001)058<1549:CTOTVO>2.0.CO;2] ]. The six cases cover a broad range of the data gathered from the datasets to validate the LES. The data used for validation comes by averaging the profiles of the field data that fall into bins corresponding to each simulation. The bins collect the meteorological tower data that is between 12LT-13LT of the corresponding month, ± 0.5σ of the corresponding surface heat flux, and ± 0.5 m / s of the hub height wind speed for the corresponding LES.

The mean velocity and temperature profiles (not shown) matched those expected for a well-mixed convective boundary layer. The potential temperature increases near the surface for the first 100-200 m. After, the potential temperature is constant up to the inversion height where there is a sharp increase in potential temperature. The increase represents the capping inversion. Finally, the potential temperature increases at a constant rate to the top boundary. The velocity profiles have a gradient for the first 100 m. Then there is a nearly constant velocity profile to the capping inversion. At the capping inversion, there is a sharp increase in velocity to the geostrophic wind condition. These match the expected mean profile of the convective boundary layer as given by Moeng [4Moeng C-H. A large-eddy-simulation model for the study of planetary boundary-layer turbulence. J Atmos Sci 1984; 41: 2052-62.
[http://dx.doi.org/10.1175/1520-0469(1984)041<2052:ALESMF>2.0.CO;2] ].

Fig. (7) compares the TKE, shear exponent, wind speed, and u^* obtained from LES to the NREL July 12LT data. The focus of the comparison is on the trends because the LES simulations are highly idealized. Fig. (7a) shows the correlation between surface heat flux and TKE at 74 m for the six LES simulations and the data (with a trendline). The LES simulations captured the TKE with an underprediction by less than 5% compared to the trendline. Fig. (7b) shows the correlation between surface heat flux and the friction velocity. LES overpredicted the friction velocities when compared to the trendline. The overprediction ranges from 8% to 16%, but the bias is relatively constant. The might be better estimated through a better determination of the friction length near the surface.

Fig. (7) (a) Correlation for surface heat flux and TKE for the positive heat flux NREL July data with trendline and LES at 74 m. (b) Correlation for surface heat flux and friction velocity for the positive heat flux NREL July data with trendline and LES at 74 m. (c) Scatter plot for surface heat flux and wind speed for the positive heat flux NREL July data and LES at 74 m. (d) Scatter for surface heat flux and shear component for the positive heat flux NREL July data and LES at 10 and 86 m.

Fig. (7c) shows the correlation between surface heat flux and wind speed at 74 m. Although they fell within the data, the estimated geostrophic wind tended to overpredict the wind speed at 74 m. The overprediction comes from the assumption of using the same values by found experimentally by Clarke and Hess [26Clarke RH, Hess GD. Geostrophic departure and the functions A and B of Rossby-number similarity theory. Boundary-Layer Meteorol 1974; 7: 267-87.
[http://dx.doi.org/10.1007/BF00240832] ] for the A and B functions values described above (more correctly would be to calculate these values for the given datasets which are outside the scope of this work). Fig. (7d) shows the correlation between surface heat flux and wind shear. The wind shear is calculated from the wind speeds between 10 m and 86 m. The over predicted wind speed creates a similar overprediction of wind shear. Again, it is noted that the wind shear falls in the range of values measured, but tends toward the higher limit of those measured values. Fig. (7) validates the methodology of using surface heat flux from experimental data for wind turbine simulations, especially for capturing the TKE near hub height.

Further validation was performed by comparing the velocity and TKE profiles for each of the simulations generated. The NREL NWTC data is collected into bins and averaged to create the mean velocity and TKE profiles. The bins collect the meteorological tower data that is between 12LT-13LT of the corresponding month, ± 0.5σ of the corresponding surface heat flux, and ± 0.5 m/s of the hub height wind speed for the corresponding LES. Fig. (8a) shows the velocity profiles for the January 12LT-1 standard deviation bin averaged data compared to the LES output. The LES output overpredicts the velocity below 50 m. Figs. (8b and 8c) shows the velocity profiles for the January mean and +1 standard deviation (respectively) bin averaged data compared to the LES output. The LES output for these two cases better predict the velocities below 50 m. The simulations capture the near uniform profiles above 50 m height that are typical for convective profiles.

Figs. (8d and 8e) show the velocity profiles for the July 12LT-1 standard deviation and mean (respectively) bin averaged data compared to the LES output. The LES output for these two cases predict the velocities below 50 m but overpredict the velocities above 50 m. Fig. (8f) shows the velocity profiles for the July 12LT +1 standard deviation bin averaged data compared to the LES output. The LES output for this case underpredicts all velocities measured from the field. Table 2 shows the error for the LES profiles compared to the binned average field data. In all simulations, the LES captures the velocity profiles well as the error falls below 5%. The small error validates the methodology for statistically selecting the boundary conditions for LES to represent day-to-day variations.

Fig. (8) Velocity profile comparisons for LES and NREL bin averaged datasets for (a) Jan -1SD, (b) Jan Mean, (c) Jan +1SD, (d) Jul -1SD, (e) Jul Mean, and (f) Jul +1SD.

Table 2
LES error for the velocity and TKE profiles.

Similarly, the TKE profiles for the bin-averaged NREL Dataset and the LES simulations are compared in Fig. (9). Fig. (9a) shows the TKE profiles for the January 12LT-1 standard deviation bin averaged data compared to the LES output. The LES underpredicts the values above 100 m. Fig. (9b) shows the TKE profiles for the January 12LT mean bin averaged data compared to the LES output. The LES underpredicts the values above 100 m and overpredicts TKE below 50 m. Fig. (9c) shows the TKE profiles for the January 12LT +1 standard deviation bin averaged data compared to the LES output. The LES overpredicts the values below 50 m but shows good agreement with all other points.

Fig. (9) Velocity profile comparisons for LES and NREL bin averaged datasets for (a) Jan -1SD, (b) Jan Mean, (c) Jan +1SD, (d) Jul -1SD, (e) Jul Mean, and (f) Jul +1SD.

Fig. (9d) shows the TKE profiles for the July 12LT -1 standard deviation bin averaged data compared to the LES output. The LES underpredicts the values above 100 m. Figs. (9e and f) show the TKE profiles for the January mean and +1 standard deviation bin averaged data compared to the LES output. The LES captures the near surface TKE well for both cases with less error above 100 m. The January 12 LT -1 standard deviation, January 12 LT Mean, and July 12 LT -1 standard deviation show small errors except for the two measurements above 100 m where the NREL data increases while the LES simulation continues to decay. The January 12LT +1 standard deviation shows a significant error at the 15 m height. The errors for the profiles are in Table 2. All errors are below 10%. The average error is 5.9% and validates that the statistically selected boundary conditions capture the TKE in the field measurements. The small error of the profiles shows that the first and second moments are obtained by using the PDF inputs for surface heat flux which further validates the methodology of creating inputs for discrete LES simulations that represent day-to-day variations in atmospheric stability.

Due to the computational cost of LES, it is common for wind energy applications to simulate a quasi-steady state of the atmospheric boundary layer with a specific atmospheric stability regime [37Churchfield MJ, Lee S, Michalakes J, Moriarty PJ. A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics. J Turbul 2012; N14.
[http://dx.doi.org/10.1080/14685248.2012.668191] ]. For example, the convective ABL simulations are often conducted by selecting a surface heat flux boundary condition that characterizes a moderately convective ABL [38Ghaisas N S, Archer C L, Xie S, Wu S, Maguire E. Evaluation of layout and atmospheric stability effects in wind farms using large-eddy simulation. Wind Energy 2017.
[http://dx.doi.org/10.1002/we.2091] -40Abkar M, Porté-Agel F. Influence of atmospheric stability on wind-turbine wakes: A large-eddy simulation study. Phys Fluids 2015; 27: 035104.
[http://dx.doi.org/10.1063/1.4913695] ]. With reducing computational cost, recent work has been conducted in simulating the ABL for the entire diurnal cycle by using the hourly values of surface heat flux obtained from field measurements [16Kumar V, Kleissl J, Meneveau C, Parlange MB. Large-eddy simulation of a diurnal cycle of the atmospheric boundary layer: Atmospheric stability and scaling issues. Water Resour Res 2006; 42: W06D09.
[http://dx.doi.org/10.1029/2005WR004651] , 18Basu S, Vinuesa J-F, Swift A. Dynamic LES modeling of a diurnal cycle. J Appl Meteorol Climatol 2008; 47: 1156-74.
[http://dx.doi.org/10.1175/2007JAMC1677.1] ]. It is becoming increasingly possible to conduct LES for the entire diurnal cycle even for wind turbine applications. LES was conducted for an idealized finite wind farm for a complete diurnal cycle [20Abkar M, Sharifi A, Porté-Agel F. Wake flow in a wind farm during a diurnal cycle. J Turbul 2016; 17: 420-41.
[http://dx.doi.org/10.1080/14685248.2015.1127379] ], and the results showed that the power deficit increased from 28% in the most convective conditions to 58% under the most stable conditions. However, it is still not computationally feasible to perform LES to obtain a monthly energy estimate required for evaluating wind potential of wind farm sites [41Celik AN. A statistical analysis of wind power density based on the Weibull and Rayleigh models at the southern region of Turkey. Renew Energy 2004; 29: 593-604.
[http://dx.doi.org/10.1016/j.renene.2003.07.002] -43García-Bustamante E, González-Rouco JF, Jiménez PA, Navarro J, Montávez JP. A comparison of methodologies for monthly wind energy estimation. Wind Energy Int J Prog Appl Wind Power Convers Technol 2009; 12: 640-59.]. In particular, variations in atmospheric stability at larger timescales (i.e., across multiple days in a month) are not always accounted for during modeling [44Barthelmie RJ, Pryor SC, Frandsen ST, et al. Quantifying the impact of wind turbine wakes on power output at offshore wind farms. J Atmos Ocean Technol 2010; 27: 1302-17.
[http://dx.doi.org/10.1175/2010JTECHA1398.1] ].

The current study proposes a novel direction for accounting the day-to-day variations in the atmospheric stability into the LES model. The variations are simulated as separate discrete LES. This methodology is similar to Monte Carlo probability theory. Monte Carlo probability uses PDFs from data to generate inputs used for discrete simulations [45James F. Monte Carlo theory and practice. Rep Prog Phys 1980; 43: 1145.
[http://dx.doi.org/10.1088/0034-4885/43/9/002] , 46Metropolis N, Ulam S. The Monte Carlo method. J Am Stat Assoc 1949; 44(247): 335-41.
[http://dx.doi.org/10.1080/01621459.1949.10483310] [PMID: 18139350] ]. The outputs of the discrete simulations help quantify the range of possible outcomes based on the variability of inputs [47Hrafnkelsson B, Oddsson GV, Unnthorsson R. A Method for estimating annual energy production using Monte Carlo Wind speed simulation. Energies 2016; 9: 286.
[http://dx.doi.org/10.3390/en9040286] ]. Although the current study does not generate the number of LES required for a complete Monte Carlo simulation [48Binder K. Introduction: Theory and “technical” aspects of Monte Carlo simulations.Monte Carlo Methods in Statistical Physics 1986; 1-45.
[http://dx.doi.org/10.1007/978-3-642-82803-4_1] ], the idea is still the same. The PDF is used to generate discrete inputs for LES (selected as the mean +1 standard deviation, and -1 standard deviation). The discrete outputs from the LES help quantify the variability in the outputs base on the variability of the inputs. The study uses high-resolution data from the NREL NWTC database to build a monthly averaged PDF of surface heat flux at each hour of the day for each month. The PDF is used to select the LES inputs at each hour, particularly the mean, +1 standard deviation, and -1 standard deviation and their corresponding initial inversion heights and geostrophic winds. The inversion height was selected based on the correlation between surface heat flux and ABL height as shown in Fig. (5). The geostrophic wind was selected using the geostrophic departure function from Clarke and Hess [26Clarke RH, Hess GD. Geostrophic departure and the functions A and B of Rossby-number similarity theory. Boundary-Layer Meteorol 1974; 7: 267-87.
[http://dx.doi.org/10.1007/BF00240832] ]. The PDF is also used to collect the field data into bins for validation. The study focuses on the PDF for July 12 LT because it had the largest mean surface heat flux and standard deviation from the mean. The study also included the January 12 LT PDF for comparisons because it had a much smaller mean surface heat flux and standard deviation.

The best fit PDF for the convective portion of the ABL was the normal distribution, which matches models of the surface heat flux for near-neutral conditions [49Katul GG. A model for sensible heat flux probability density function for near-neutral and slightly-stable atmospheric flows. Boundary-Layer Meteorol 1994; 71: 1-20.
[http://dx.doi.org/10.1007/BF00709217] ]. The mean of the PDF is the expected value of the function over the entire month at each hour. The standard deviation represents the variation that occurs at each hour of the month. The difference between the +1 standard deviation and -1 standard deviation represents 68% of the data. The field data was collected into bins that correspond to each simulation. The bins collect the meteorological tower data that is between 12 -13LT of the corresponding month, ± 0.5 σ of the corresponding surface heat flux, and ± 0.5 m / s of the hub height wind speed for the corresponding LES. The binned metrological data were averaged to generate vertical profiles of velocity and TKE.

LES was conducted to simulate the case of January 12LT and July 12LT by using the mean value of PDFs of surface flux as discrete boundary conditions. For all the simulations, the initial inversion height was selected using the correlation between surface heat flux and boundary height measured in the field data. The geostrophic wind was estimated using u^* with the geostrophic departure described by Clarke and Hess [26Clarke RH, Hess GD. Geostrophic departure and the functions A and B of Rossby-number similarity theory. Boundary-Layer Meteorol 1974; 7: 267-87.
[http://dx.doi.org/10.1007/BF00240832] ]. The mean July simulation had a higher hub height wind speed (5.0vs. 4.2 m/s), and TKE (1.4vs. 0.7 m²/s²) than the January mean simulation. Two additional LES were conducted by using -1 standard deviation, and + 1 standard deviation from mean values of PDF of surface flux as inputs into LES. The July +1 standard deviation simulation had a hub height TKE of 1.7 m²/s² while the -1 standard deviation simulation was 0.9 m²/s². The increased surface heat flux increased the hub height TKE by 2 times. The difference in the January +1 standard deviation and -1 standard deviation was 0.35 m²/s² and 0.95 m²/s². The day-to-day variations in surface heat flux had significant impacts on the outputs of LES. These wide variations in TKE flux have considerable impacts on wind turbine performance [22Wu Y-T, Porté-Agel F. Atmospheric turbulence effects on wind-turbine wakes: An LES study. Energies 2012; 5: 5340-62.
[http://dx.doi.org/10.3390/en5125340] , 24Clifton A, Schreck S, Scott G, Kelley N, Lundquist JK. Turbine inflow characterization at the national wind technology center. J Sol Energy Eng 2013; 135: 031017-031017-11.
[http://dx.doi.org/10.1115/1.4024068] , 36Iungo GV, Porté-Agel F. Volumetric lidar scanning of wind turbine wakes under convective and neutral atmospheric stability regimes. J Atmos Ocean Technol 2014; 31: 2035-48.
[http://dx.doi.org/10.1175/JTECH-D-13-00252.1] ].

The outputs of the six LES captured the same correlation of surface heat flux with TKE, wind speed, wind shear, and friction velocity at hub height as the field measurements. The LES correlation between surface heat flux and TKE showed that for every increase in 0.1 kms^-1 the TKE increased by 0.42 m²/s². The LES correlations showed an increase in friction velocity and wind speed of 0.095 m/s and 0.74 m/s for every increase in 0.1 kms^-1. The correlations show the significant impact to LES outputs due to the variation in surface heat flux. Using the methodology described the discrete LES were able to capture these variations using discrete LES.

The binned average velocity and TKE profiles validated the LES. The LES profiles matched the velocity profiles with 4% error and the TKE profiles to 6% error on average. The results are promising and demonstrate the methodology approach is robust enough to capture the ABL turbulence over multiple diurnal cycles. The current work is an important direction for futuristic LES studies of wind turbines in ABL based on field measurements. The significance of the results is two folds (a) Day-to-day variability of surface heat flux is well-represented in TKE and ABL depth, hence surface heat flux is a boundary condition that dictates the efficacy of LES in capturing the day-to-day variation of ABL turbulence, and (b) The methodology allows the outputs from LES to represent mathematically significant variations that occur in the atmosphere. Applications of LES, such as wind turbine simulations, can use the methodology to quantify the variability in turbine performance based on the day-to-day variations of surface heat flux. These predictions add further fidelity to wind turbine LES. The study serves as a validation incorporate field measurements averaged over multiple diurnal cycles into LES applications.