Javier Sanz Rodrigo, National Renewable Energy Centre (CENER), Sarriguren, Spain, jsrodrigo@cener.com
July 2017
This benchmark report analyzes the simulations submitted to the GABLS3 benchmark for wind energy, organized within the NEWA project and the IEA Task 31 Wakebench. The results of this benchmark have been published in:
Sanz Rodrigo J., Allaerts D., Ávila M., Barçons J., Cavar D., Chávez Arroyo R.A., Churchfield M., Draper M., Kosovic B., Lundquist J.K., Meyers J., Muñoz Esparza D., Palma J.M.L.M., Tomazewski J.M., Troldborg N., van der Laan M.P. and Veiga Rodrigues C. (2017) Results of the GABLS3 diurnal cycle benchmark for wind energy applications. Journal of Physics: Conference Series, 854: 012037, https://doi.org/10.1088/1742-6596/854/1/012037
Background information and benchmark set-up can be found in: http://windbench.net/gabls-3
The following simulations participate in the benchmark. Notice that WRF-YSU (ref) is the reference mesoscale simulation that was used to generate the mesoscale tendencies that are used as forcings for the microscale models. Additional WRF simulations have been run to test the sensitivity of mesoscale simulations to input forcing and planetary boundary layer (PBL) scheme. The ensemble mean of these WRF simulations is also used in the model intercomparison.
Table 1. Summary of model simulations. Monin Obukhov similarity theory (Monin-Obukhov (MOST) surface boundary conditions use either heat flux ($H$), 2-m ($T_{2}$) or skin temperature ($T_{SK}$) from WRF.
| Name | Input | Turbulence | z-Levels | Surface B.C. |
|---|---|---|---|---|
| WRF-YSU (ref) | ERA-Interim | YSU | 46 | Noah |
| WRF | ERA-Interim, GFS | MJY, MYNN, QNSE, TEMF, YSU | 46 | Noah |
| WRF-YSU_LES | ERA-Interim | LES-TKE | 101 | Noah |
| WRF-VentosM_ke | ERA-Interim | YSU/$k-\epsilon$ | 70 | MOST, H |
| CFDWind1D_ke | WRF (ref) | $k-\epsilon$ | 301 | MOST, $T_2$ |
| Alya-CFDWind1D_ke | WRF (ref) | $k-\epsilon$ | 500 | MOST, $T_2$ |
| Ellipsys1D_ke | WRF (ref) | $k-\epsilon$ | 512 | MOST, $T_2$ |
| Ellipsys3D_ke | WRF (ref) | $k-\epsilon$ | 192 | MOST, $T_{SK}$ |
| Ellipsys3D_LES | WRF (ref) | Smagorinsky | 128 | MOST, $T_{SK}$ |
| SP-Wind_LES | WRF (ref) | LES-TKE | 500 | MOST, $T_2$ |
In [11]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy import spatial
import datetime
import netCDF4
import pandas as pd
from scipy import interpolate
import scipy.integrate as integrate
import statsmodels.api as sm
import os
# Constants
K = 0.41 # von Karman constant
g = 9.81 # [m s-2]
R_air = 287.058 # Specific gas constant for dry air [J kg-1 K-1]
Cp_air = 1005 # Specific heat of air [J kg-1 K-1]
P0 = 100000. # Reference pressure [Pa]
# Site ID
siteID = 'GABLS3'
lat_s = 51.971 # degrees N
lon_s = 4.927 # degrees E
# Evaluation Period and reference Rotor chararcteristics
datefrom = datetime.datetime(2006, 7, 1, 12, 0, 0) # Origin of evaluation period
dateto = datetime.datetime(2006, 7, 2, 12, 0, 0) # End of evaluation period
Hhub = 120. # hub-height
Drot = 160. # rotor diameter
ts = 10 # sampling frequency to evaluate [min]
# Mesoscale tendencies averaging settings
tav = 60.0 # Time averaging time used in simulations [min]
Lav = 9000.0 # Spatial averaging [m]
ws = int(tav/ts) # rolling window size
In [12]:
dirsim = './data'
filesim = [dirsim + '/GABLS3_tendencies_d02_YSU_w60_L9000.nc',
dirsim + '/WRF-LES_d8.nc',
dirsim + '/VentosM.nc',
dirsim + '/CFDWindSCM_ke.nc',
dirsim + '/Alya-CFDWind1D.nc',
dirsim + '/Ellipsys1D_TskWRF.nc',
dirsim + '/Ellipsys3D_TskWRF.nc',
dirsim + '/Ellipsys3D_LES_TskWRF.nc',
dirsim + '/SP-Wind.nc',
dirsim + '/WRF-YSU_GFS.nc',
dirsim + '/WRF-YSU_ERA.nc',
dirsim + '/WRF-TEMF_GFS.nc',
dirsim + '/WRF-TEMF_ERA.nc',
dirsim + '/WRF-QNSE_GFS.nc',
dirsim + '/WRF-QNSE_ERA.nc',
dirsim + '/WRF-MYNN_GFS.nc',
dirsim + '/WRF-MYNN_ERA.nc',
dirsim + '/WRF-MJY_GFS.nc',
dirsim + '/WRF-MJY_ERA.nc',
]
simID = ['WRF-YSU (ref)','WRF-YSU_LES','VentosM',
'CFDWind1D_ke','Alya-Wind1D_ke','Ellipsys1D_ke',
'Ellipsys3D_ke',
'Ellipsys3D_LES','SP-Wind_LES',
'WRF-YSU_GFS','WRF-YSU_ERA',
'WRF-TEMF_GFS','WRF-TEMF_ERA',
'WRF-QNSE_GFS','WRF-QNSE_ERA',
'WRF-MYNN_GFS','WRF-MYNN_ERA',
'WRF-MJY_GFS','WRF-MJY_ERA',
]
simtype = ['meso','micro','micro',
'micro','micro','micro',
'micro',
'micro','micro',
'meso','meso',
'meso','meso',
'meso','meso',
'meso','meso',
'meso','meso',
]
Nsim = len(filesim)
# Create a list with the simulation datasets
t = []; U = []; V = []; Th = []; z = []; S = []; WD = []; us = []
wt = []; T2 = []; TKE = []
for isim in range(0,Nsim):
f = netCDF4.Dataset(filesim[isim], 'r')
times = f.variables['time'][:]
idates = np.where(np.logical_and(times >= mdates.date2num(datefrom), times < mdates.date2num(dateto)))[0]
t.append(pd.DataFrame(f.variables['time'][idates], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean())
U.append(pd.DataFrame(f.variables['U'][idates,:], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
V.append(pd.DataFrame(f.variables['V'][idates,:], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
Th.append(pd.DataFrame(f.variables['Th'][idates,:], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
z.append(f.variables['z'][:])
S.append((U[isim]**2 + V[isim]**2)**0.5)
WD.append(180 + np.arctan2(U[isim],V[isim])*180/np.pi)
us.append(pd.DataFrame(f.variables['ust'][idates], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
wt.append(pd.DataFrame(f.variables['wt'][idates], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
T2.append(pd.DataFrame(f.variables['T2'][idates], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
if simtype[isim] == 'micro':
TKE.append(pd.DataFrame(f.variables['TKE'][idates,:], index = mdates.num2date(f.variables['time'][idates])).resample(str(ts)+'Min').mean().rolling(window = ws).mean())
else:
TKE.append(pd.DataFrame(data=None, columns=U[isim].columns,index=U[isim].index))
f.close()
# Create ensemble from WRF simulations
simID.append('WRF-Ensemble')
simtype.append('meso')
Nsim = Nsim + 1
U.append(pd.DataFrame(data=np.mean(np.dstack(U[9:18]), axis = 2), columns=U[18].columns,index=U[18].index))
V.append(pd.DataFrame(data=np.mean(np.dstack(V[9:18]), axis = 2), columns=U[18].columns,index=U[18].index))
S.append((U[len(U)-1]**2 + V[len(U)-1]**2)**0.5)
WD.append(180 + np.arctan2(U[len(U)-1],V[len(U)-1])*180/np.pi)
Th.append(pd.DataFrame(data=np.mean(np.dstack(Th[9:18]), axis = 2), columns=U[18].columns,index=U[18].index))
TKE.append(pd.DataFrame(data=np.mean(np.dstack(TKE[9:18]), axis = 2), columns=U[18].columns,index=U[18].index))
us.append(pd.DataFrame(data=np.mean(np.dstack(us[9:18]), axis = 1), index=U[18].index))
wt.append(pd.DataFrame(data=np.mean(np.dstack(wt[9:18]), axis = 1), index=U[18].index))
T2.append(pd.DataFrame(data=np.mean(np.dstack(T2[9:18]), axis = 1), index=U[18].index))
t.append(t[18])
z.append(z[18])
In [13]:
# Note that the file 'gabls3_scm_cabauw_obs_v33.nc' can be obtained from the KNMI GABLS3 website
# http://projects.knmi.nl/gabls/gabls3_scm_cabauw_obs_v33.nc
# Alternatively, it is also provided in the Windbench/GABLS3 input dataset
dirobs = './data'
fileobs = dirobs + '/gabls3_scm_cabauw_obs_v33.nc'
nodata = -9999.0 # missing data flag
date0 = datetime.datetime(2006,7,1,0,0,0) # origin of time_obs
f1 = netCDF4.Dataset(fileobs, 'r')
dates_obs = f1.variables['date'][:]
time_obs = f1.variables['time'][:] # 'hours since 2006-07-01 00:00:00 0:00'
date_obs = []
for i in range (0,len(time_obs)):
date_obs.append(date0 + datetime.timedelta(seconds = np.int(time_obs[i]*3600.0)))
ifrom_obs=0
for j in range(0,len(date_obs)):
if date_obs[j] <= datefrom:
ifrom_obs = j
ito_obs=0
for j in range(0,len(date_obs)):
if date_obs[j] <= dateto:
ito_obs = j+1
date_obs = mdates.date2num(date_obs[ifrom_obs:ito_obs])
time_obs = 24.0*(date_obs - date_obs[0])
# Surface variables
ustar_obs = pd.DataFrame(np.ma.array(f1.variables['ustar'][ifrom_obs:ito_obs], mask=(f1.variables['ustar'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
U10_obs = pd.DataFrame(np.ma.array(f1.variables['u10m'][ifrom_obs:ito_obs], mask=(f1.variables['u10m'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
V10_obs = pd.DataFrame(np.ma.array(f1.variables['v10m'][ifrom_obs:ito_obs], mask=(f1.variables['v10m'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
S10_obs = pd.DataFrame(np.ma.array(f1.variables['f10m'][ifrom_obs:ito_obs], mask=(f1.variables['f10m'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
WD10_obs = pd.DataFrame(np.ma.array(f1.variables['d10m'][ifrom_obs:ito_obs], mask=(f1.variables['d10m'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
T2_obs = pd.DataFrame(np.ma.array(f1.variables['t2m'][ifrom_obs:ito_obs], mask=(f1.variables['t2m'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
HFX_obs = pd.DataFrame(np.ma.array(f1.variables['shf'][ifrom_obs:ito_obs], mask=(f1.variables['shf'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
Ts_obs = pd.DataFrame(np.ma.array(f1.variables['ts'][ifrom_obs:ito_obs], mask=(f1.variables['ts'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
zs_obs = f1.variables['zs'][ifrom_obs:ito_obs][0,:] # soil temperature heights [m]
Tsk_obs = pd.DataFrame(np.ma.array(f1.variables['tsk'][ifrom_obs:ito_obs], mask=(f1.variables['tsk'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
Tsk1_obs = pd.DataFrame(np.ma.array(f1.variables['tsk1'][ifrom_obs:ito_obs], mask=(f1.variables['tsk1'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
# Vertical profiles
zf_obs = f1.variables['zf'][ifrom_obs:ito_obs][0,:] # velocity profile heights [m]
zT_obs = f1.variables['zt'][ifrom_obs:ito_obs][0,:] # temperature profile heights [m]
U_obs = pd.DataFrame(np.ma.array(f1.variables['u'][ifrom_obs:ito_obs], mask=(f1.variables['u'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
V_obs = pd.DataFrame(np.ma.array(f1.variables['v'][ifrom_obs:ito_obs], mask=(f1.variables['v'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
S_obs = pd.DataFrame(np.ma.array(f1.variables['f'][ifrom_obs:ito_obs], mask=(f1.variables['f'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs)) # horizontal wind speed [deg]
WD_obs = pd.DataFrame(np.ma.array(f1.variables['d'][ifrom_obs:ito_obs], mask=(f1.variables['d'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs)) # wind direction [deg]
Th_obs = pd.DataFrame(np.ma.array(f1.variables['th'][ifrom_obs:ito_obs], mask=(f1.variables['th'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs)) # potential temperature [K]
q_obs = pd.DataFrame(np.ma.array(f1.variables['q'][ifrom_obs:ito_obs], mask=(f1.variables['q'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
# fluxes
zh_obs = f1.variables['zh'][ifrom_obs:ito_obs][0,:] # fluxes heights
wt_obs = pd.DataFrame(np.ma.array(f1.variables['wt'][ifrom_obs:ito_obs], mask=(f1.variables['wt'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
uw_obs = pd.DataFrame(np.ma.array(f1.variables['uw'][ifrom_obs:ito_obs], mask=(f1.variables['uw'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
vw_obs = pd.DataFrame(np.ma.array(f1.variables['vw'][ifrom_obs:ito_obs], mask=(f1.variables['vw'][ifrom_obs:ito_obs] == nodata)), index = mdates.num2date(date_obs))
us_obs = (uw_obs**2 + vw_obs**2)**0.25
# Obukhov length
T0_obs = np.tile(Th_obs[36].values,(5,1)).T
L_obs = -us_obs**3/(K*(g/T0_obs)*wt_obs)
# Resample and apply tav to observations
ws = int(tav/ts) # rolling window size
ustar_obs_w = ustar_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
HFX_obs_w = HFX_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
U10_obs_w = U10_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
V10_obs_w = V10_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
S10_obs_w = S10_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
WD10_obs_w = WD10_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
T2_obs_w = T2_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
Ts_obs_w = Ts_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
U_obs_w = U_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
V_obs_w = V_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
S_obs_w = S_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
WD_obs_w = WD_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
Th_obs_w = Th_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
q_obs_w = q_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
wt_obs_w = wt_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
uw_obs_w = uw_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
vw_obs_w = vw_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
us_obs_w = us_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
L_obs_w = L_obs.resample(str(ts)+'Min').mean().bfill().rolling(window = ws).mean()
TKE_obs_w = pd.DataFrame(data=None, columns=U_obs_w.columns,index=U_obs_w.index)
In [14]:
# List simulations for index reference
for isim in range(0,len(simID)):
print isim, simID[isim]
In [15]:
# Select the simulations you want to plot
plotsim = np.array([0,19,1,2,3,4,5,6,7,8])
# Plot settings
zlim = 2000 # height limit [m]
Zcmap = plt.get_cmap('jet') # colormap
figname = siteID+'_tzfields.png' # output filename
Zlevels = np.array([np.linspace(0,16,17), # range for S
np.linspace(30,210,13), # range for WD
np.linspace(285,305,11), # range for Th
np.linspace(0,5,11)]) # range for TKE
X = []; Y = []; Z = []; taxis = []; zaxis = []; Ztitle = []
Z = [S_obs_w.T, WD_obs_w.T, Th_obs_w.T, TKE_obs_w.T]
[Xf,Yf] = np.meshgrid(date_obs,zf_obs)
[XT,YT] = np.meshgrid(date_obs,zT_obs)
X = [Xf, Xf, XT, Xf]
Y = [Yf, Yf, YT, Yf]
taxis = [date_obs, date_obs, date_obs, date_obs]
zaxis = [zf_obs, zf_obs, zT_obs, zf_obs]
Ztitle = ['Observations']
for ip in range(0,len(plotsim)):
isim = plotsim[ip]
Z.append(S[isim].T)
[X0,Y0] = np.meshgrid(t[isim],z[isim])
taxis.append(t[isim]); zaxis.append(z[isim]); X.append(X0); Y.append(Y0)
Z.append(WD[isim].T)
[X0,Y0] = np.meshgrid(t[isim],z[isim])
taxis.append(t[isim]); zaxis.append(z[isim]); X.append(X0); Y.append(Y0)
Z.append(Th[isim].T)
[X0,Y0] = np.meshgrid(t[isim],z[isim])
taxis.append(t[isim]); zaxis.append(z[isim]); X.append(X0); Y.append(Y0)
Z.append(TKE[isim].T)
[X0,Y0] = np.meshgrid(t[isim],z[isim])
taxis.append(t[isim]); zaxis.append(z[isim]); X.append(X0); Y.append(Y0)
Ztitle.append(simID[isim])
Zlabel = ('$S$ [$m s^{-1}$]','$WD$ ['+u'\N{DEGREE SIGN}'+']','$\Theta$ [$K$]','$TKE$ [$m^{2} s^{-2}$]')
taxis_label = 'UTC time in hours since ' + datefrom.strftime('%Y-%m-%d %H:%M') + ', $L_{avg}$ = ' + "%.0f"%(Lav/1000) + ' km, $t_{avg}$ = ' + "%.0f"%(tav) + ' min'
zaxis_label = '$z$ [$m$]'
hoursFmt = mdates.DateFormatter('%H')
Zlevels = np.tile(Zlevels,len(Ztitle))
hoursFmt = mdates.DateFormatter('%H')
tticks = [datefrom + datetime.timedelta(hours = 3*x) for x in range(0, 9)]
fig, ax = plt.subplots(len(Ztitle), 4, sharex='col', sharey='row', figsize=(11,17), dpi=300)
xrotor = np.array([mdates.date2num(datefrom), mdates.date2num(dateto)])
yrotor1 = np.array([Hhub - 0.5*Drot, Hhub - 0.5*Drot])
yrotor2 = np.array([Hhub + 0.5*Drot, Hhub + 0.5*Drot])
cbleft = np.array([0.125,0.33,0.53,0.73])
for iax in range (0,len(Ztitle)*4):
ix,iy = np.unravel_index(iax,(len(Ztitle),4))
CF0 = ax[ix,iy].contourf(X[iax],Y[iax],Z[iax], Zlevels[iax], cmap=Zcmap)
ax[ix,iy].plot(xrotor,yrotor1,'--k')
ax[ix,iy].plot(xrotor,yrotor2,'--k')
ax[ix,iy].set_ylim([10, zlim])
ax[ix,iy].set_yscale('log')
ax[ix,iy].xaxis.set_major_formatter(hoursFmt)
ax[ix,iy].set_xticks(tticks)
if iy == 0:
ax[ix,iy].set_ylabel(zaxis_label)
if iy == 3:
ax[ix,iy].yaxis.set_label_position("right")
ax[ix,iy].set_ylabel(Ztitle[ix], fontsize = 9, fontweight='bold')
if ix == len(Ztitle)-1:
fig.subplots_adjust(top=0.87)
cbar_ax = fig.add_axes([cbleft[iy],0.885, 0.165, 0.012])
cbar = fig.colorbar(CF0, cax=cbar_ax, orientation='horizontal')
cbar.ax.set_xlabel(Zlabel[iy],labelpad = -42, x = 0.5)
ax[len(Ztitle)-1,1].set_xlabel(taxis_label,x=1.2);
plt.setp([a.get_xticklabels() for a in ax[0, :]], visible=False)
plt.setp([a.get_yticklabels() for a in ax[:, 1]], visible=False)
#fig.savefig(figname, bbox_inches='tight', dpi = 300)
plt.show()
from IPython.display import Markdown, display
figcaption = ("**Fig 5. Time-height contour of horizontal wind speed $S$, direction $WD$, "
"potential temperature $\Theta$ and turbulent kinetic energy $TKE$. Dotted lines "
"denote a rotor diameter of 160 m at 120 m hub-height.**")
display(Markdown(figcaption))
In [16]:
# Specify the datetime at which you want to plot the vertical profiles
t0 = datetime.datetime(2006,7,2,0,0,0)
# Plot settings
linespec = ['k-','b.-','r.-','c.-','m-','g-','y-','c-','b-.','r--']
lwidth = np.array([2,1,1,1,1,1,1,1,2,2])
Nm = 3 # marker every
Z_obs = (S_obs_w.loc[t0].values, WD_obs_w.loc[t0].values,
Th_obs_w.loc[t0].values)
z_obs = (zf_obs, zf_obs, zT_obs)
Z_sim = []; z_sim = []
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
Z_sim.append((S[isim].loc[t0].values, WD[isim].loc[t0].values,
Th[isim].loc[t0].values))
z_sim.append((z[isim], z[isim], z[isim]))
figname = siteID+'_'+t0.strftime('%Y-%m-%d_%H')+'_profiles.png'
zlim = 1000
fig,ax = plt.subplots(1, 3, sharey='row', figsize=(8,6))
for iax in range (0,3):
ax[iax].plot(Z_obs[iax], z_obs[iax], 'ok', color='grey', label = 'obs ')
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
if (isim == 0 and iax == 2):
aa = []
else:
ax[iax].plot(Z_sim[iplot][iax][1:], z_sim[iplot][iax][1:],
linespec[iplot], linewidth = lwidth[iplot], label = simID[isim], markevery=Nm)
#ax[iax].set_yscale('log')
xlim = ax[iax].get_xlim()
ax[iax].plot(xlim,yrotor1,'--k')
ax[iax].plot(xlim,yrotor2,'--k')
ax[iax].set_ylim([1., zlim])
ax[iax].set_yticks(np.linspace(0,zlim,5))
ax[iax].grid(which='major',color='k',linestyle=':')
#ax[0].set_xlim([0, 16])
ax[1].set_xlim([60, 180])
ax[2].set_xlim([285, 300])
ax[0].set_xlabel('$S$ [$m s^{-1}$]')
ax[1].set_xlabel('$WD$ ['+u'\N{DEGREE SIGN}'+']')
ax[2].set_xlabel(r'$\Theta$ [$K$]')
ax[0].set_ylabel('$z$ [$m$]')
ax[0].set_title('$t_{0}$ = '+ t0.strftime('%Y-%m-%d %H:%M:%S'))
ax[1].legend(prop={'size':10},bbox_to_anchor=(-0.5, 1), loc='upper left')
plt.setp([a.get_yticklabels() for a in ax[1:]], visible=False)
plt.tight_layout()
#fig.savefig(figname, bbox_inches='tight', format='png', dpi=300)
plt.show()
figcaption = ("**Fig 7. Vertical profiles of wind speed, wind direction "
"and potential temperature at UTC **"+t0.strftime('%Y-%m-%d %H:%M:%S'))
display(Markdown(figcaption))
The performance of the models is based on quantities of interest relevant for wind energy applications. These quantities are evaluated across a reference rotor span of 160 m, between 40 and 200 m, characteristic of an 8-MW large wind turbine. Besides hub-height wind speed $S_{hub}$ and direction $WD_{hub}$, it is relevant to consider the rotor-equivalent wind speed $REWS$, the turbulence intensity (not evaluated here), the wind speed shear $\alpha$, and the wind direction shear or veer $\psi$.
The REWS is especially suitable to account for wind shear in wind turbine power performance tests [14]. The REWS is the wind speed corresponding to the kinetic energy flux through the swept rotor area, when accounting for the vertical shear:
$$ REWS = \left[\frac{1}{A}\sum_{i}(A_iS_i^3\cos\beta_i)\right]^{1/3} $$where $A$ is the rotor area and $A_i$ are the horizontal segments that separate vertical measurement points of horizontal wind speed $S_i$ across the rotor plane. The $REWS$ is here weighted by the cosine of the angle $\beta_i$ of the wind direction $WD_i$ with respect to the hub-height wind direction to account for the effect of wind veer [15].
Wind shear is defined by fitting a power-law curve across the rotor wind speed points $S_i$:
$$ S_i = S_{hub}\left(\frac{z_i}{z_{hub}}\right)^{\alpha} $$Similarly, wind veer is defined as the slope $\psi$ of the linear fit of the wind direction difference:
$$ \beta_i = \psi(z_i-z_{hub}) $$To evaluate simulations and measurements consistently, these quantities are obtained after resampling, by linear interpolation, velocity and wind direction vertical profiles at 10 points across the rotor area and then computing the REWS and the shear functional fits. While these fitting functions are commonly used in wind energy, their suitability in LLJ conditions is questionable. The regression coefficient from the fitting can be used to determine this suitability.
A rolling average with a window size of one hour is applied to simulation and observational data to remove the impact of high frequency fluctuations in the analysis.
Validation results are quantified in terms of the mean absolute error ($MAE$):
$$ MAE = \frac{1}{N}\sum_{i=1}^{N}\left|\chi_{pred}-\chi_{obs}\right| $$where $\chi$ is any of the above mentioned quantities of interest, predicted (pred) or observed (obs), and $N$ is the number of samples evaluated in the time series.
In [17]:
# Define rotor-based quantities of interest
zrot = np.linspace(Hhub - 0.5*Drot, Hhub + 0.5*Drot, 1 + Drot/10, endpoint=True)
def rotor(z,Sz,WDz):
# Rotor QoIs [m s-1]
# z: heights [m] where the velocity and wind direction are known spanning the rotor diameter
# Sz, WDz: Wind speed [m s-1] and direction [deg from N] at z levels [tdim,zdim]
# Returns:
# REWS: rotor equivalent wind speed [m s-1]
# alpha: wind shear, power-law exponent from linear fit lg(U/Uhub) = alpha*log(z/zhub)
# alpha_R2: R-squared from least squares fit to linear function
# veer: wind veer, slope of linear function beta = WDz - WDhub = veer*(z - zhub)
# veer_R2: R-squared from least squares fit to linear function
tdim = Sz.shape[0]
zdim = Sz.shape[1]
Rrot = 0.5*(z[-1]-z[0])
Hhub = 0.5*(z[-1]+z[0])
ihub = int(0.5*len(zrot))
Arotor = np.pi*(Rrot)**2
Uz = -Sz*np.sin(np.pi*WDz/180.)
Vz = -Sz*np.cos(np.pi*WDz/180.)
Shub = Sz[:,ihub]
WDhub = WDz[:,ihub]
def cz(x,R,H):
return 2.*(R**2 - (x-H)**2)**0.5
sumA = np.zeros((Sz.shape[0]))
veer = np.zeros((Sz.shape[0]))
for i in range(0,zdim-1):
Ai = integrate.quad(cz, z[i], z[i+1], args = (Rrot,Hhub))
Si = 0.5*(Sz[:,i+1]+Sz[:,i])
Ui = 0.5*(Uz[:,i+1]+Uz[:,i])
Vi = 0.5*(Vz[:,i+1]+Vz[:,i])
WDi = 180. + np.arctan2(Ui,Vi)*180./np.pi
betai = WDi - WDhub
sumA = sumA + Ai[0]*(Si*np.cos(np.pi*betai/180.))**3
REWS = (sumA/Arotor)**(1./3.)
alpha = np.zeros(tdim); alpha_stderr = np.zeros(tdim); alpha_R2 = np.zeros(tdim)
veer = np.zeros(tdim); veer_stderr = np.zeros(tdim); veer_R2 = np.zeros(tdim)
for it in range(0,tdim):
regmodel = sm.OLS(np.log(Sz[it,:]/Shub[it]), np.log(z/Hhub))
results = regmodel.fit()
alpha[it] = results.params[0]
alpha_stderr[it] = results.bse[0]
alpha_R2[it] = results.rsquared
regmodel = sm.OLS(WDz[it,:] - WDhub[it], z - Hhub)
results = regmodel.fit()
veer[it] = results.params[0]
veer_stderr[it] = results.bse[0]
veer_R2[it] = results.rsquared
return REWS, Shub, WDhub, alpha, alpha_R2, veer, veer_R2
Srews_obs = interpolate.interp2d(zf_obs[30:34],date_obs,S_obs_w.values[:,30:34])(zrot,date_obs)
Urews_obs = interpolate.interp2d(zf_obs[30:34],date_obs,U_obs_w.values[:,30:34])(zrot,date_obs)
Vrews_obs = interpolate.interp2d(zf_obs[30:34],date_obs,V_obs_w.values[:,30:34])(zrot,date_obs)
WDrews_obs = 180. + np.arctan2(Urews_obs,Vrews_obs)*180./np.pi
REWS_obs, Shub_obs, WDhub_obs, alpha_obs, alpha_R2_obs, veer_obs, veer_R2_obs = rotor(zrot,Srews_obs,WDrews_obs)
REWS = []; Shub = []; WDhub = []
alpha = []; alpha_R2 = []; veer = []; veer_R2 = []
for isim in range(0,Nsim):
Srews_sim = interpolate.interp2d(z[isim],np.ravel(t[isim].values),S[isim].values)(zrot,np.ravel(t[isim].values))
Urews_sim = interpolate.interp2d(z[isim],np.ravel(t[isim].values),U[isim].values)(zrot,np.ravel(t[isim].values))
Vrews_sim = interpolate.interp2d(z[isim],np.ravel(t[isim].values),V[isim].values)(zrot,np.ravel(t[isim].values))
WDrews_sim = 180. + np.arctan2(Urews_sim,Vrews_sim)*180./np.pi
REWS0, Shub0, WDhub0, alpha0, alpha_R20, veer0, veer_R20 = rotor(zrot,Srews_sim,WDrews_sim)
REWS.append(REWS0)
Shub.append(Shub0)
WDhub.append(WDhub0)
alpha.append(alpha0)
alpha_R2.append(alpha_R20)
veer.append(veer0)
veer_R2.append(veer_R20)
In [18]:
# plot QoIs
figname = siteID+'_QoIs_rotor.png'
taxis_label = 'UTC time in hours since ' + datefrom.strftime('%Y-%m-%d %H:%M') + ', $D$ = ' + "%.0f"%(Drot) + ' m, $H$ = ' + "%.0f"%(Hhub) + ' m'
fig, ax = plt.subplots(2, 2, sharex='col', figsize=(9,6), dpi=300)
ax[0][0].plot(date_obs[::3], REWS_obs[::3], 'ok', color='grey', label = 'obs ')
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
ax[0][0].plot(t[isim], REWS[isim], linespec[iplot],linewidth = lwidth[iplot],
label = simID[isim], markevery=Nm)
ax[0][0].set_ylabel(r'$REWS$ [$m s^{-1}$]')
ax[0][0].grid(which='major',color='grey',linestyle=':')
ax[0][0].set_ylim([2., 14.])
ax[1][0].plot(date_obs[::3], WDhub_obs[::3], 'ok', color='grey', label = 'obs ')
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
ax[1][0].plot(t[isim], WDhub[isim], linespec[iplot],linewidth = lwidth[iplot],
label = simID[isim], markevery=Nm)
ax[1][0].set_ylabel(r'$WD_{hub}$ ['+u'\N{DEGREE SIGN}'+']')
ax[1][0].grid(which='major',color='grey',linestyle=':')
ax[1][0].set_ylim([60., 180.])
ax[1][0].xaxis.set_major_formatter(hoursFmt)
ax[1][0].set_xticks(tticks)
ax[1][0].set_xlim([mdates.date2num(datefrom), mdates.date2num(dateto)])
ax[1][0].set_xlabel(taxis_label)
ax[1][0].xaxis.set_label_coords(1.2, -0.1)
ax[0][1].plot(date_obs[::3], alpha_obs[::3], 'ok', color='grey', label = 'obs ')
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
ax[0][1].plot(t[isim], alpha[isim], linespec[iplot],linewidth = lwidth[iplot],
label = simID[isim], markevery=Nm)
ax[0][1].set_ylabel(r'Wind Shear $\alpha$')
ax[0][1].grid(which='major',color='grey',linestyle=':')
ax[0][1].set_ylim([-0.2, 1.0])
ax[1][1].plot(date_obs[::3], veer_obs[::3], 'ok', color='grey', label = 'obs ')
for iplot in range(0,len(plotsim)):
isim = plotsim[iplot]
ax[1][1].plot(t[isim], veer[isim], linespec[iplot],linewidth = lwidth[iplot],
label = simID[isim], markevery=Nm)
ax[1][1].set_ylabel(r'Wind Veer $\psi$')
ax[1][1].grid(which='major',color='grey',linestyle=':')
ax[1][1].set_ylim([-0.1, 0.6])
ax[1][1].xaxis.set_major_formatter(hoursFmt)
ax[1][1].set_xticks(tticks)
ax[1][1].set_xlim([mdates.date2num(datefrom), mdates.date2num(dateto)])
ax[0][1].legend(prop={'size':8}, bbox_to_anchor=(1.47, 1), ncol=1, loc='upper right')
fig.subplots_adjust(wspace=.25)
plt.show()
#fig.savefig(figname, bbox_inches='tight', dpi = 300)
figcaption = ("**Fig 6. Time series of rotor rotor-based quantities of interest used for validation**")
display(Markdown(figcaption))
In [28]:
from IPython.display import display
pd.options.display.float_format = '{:,.2f}'.format # format for output data in tables
# Compute error metics
def mae(xtrue, xpred, norm):
# Normalized Mean Absolute Error
# xtrue: series with true values
# xpred: series with predicted values
nonan = ~np.isnan(xtrue)
if norm == 'True':
xtrue0 = xtrue[nonan]/np.mean(xtrue[nonan])
xpred0 = xpred[nonan]/np.mean(xtrue[nonan])
else:
xtrue0 = xtrue[nonan]
xpred0 = xpred[nonan]
N = len(xtrue0)
abserr = np.abs(xpred0 - xtrue0)
MAE = np.sum(abserr)/N
return MAE
# Errors with respecto to observations
Metrics = pd.DataFrame(np.zeros((Nsim, 5)),
columns = ['REWS','Shub','WDhub','alpha','veer'],
index = simID)
Metrics.loc['units',:] = np.array(['m s-1','m s-1','deg','-','-'])
for isim in range(0,Nsim):
Metrics.loc[(simID[isim]),'REWS'] = mae(REWS_obs,REWS[isim],'False')
Metrics.loc[(simID[isim]),'Shub'] = mae(Shub_obs,Shub[isim],'False')
Metrics.loc[(simID[isim]),'WDhub'] = mae(WDhub_obs,WDhub[isim],'False')
Metrics.loc[(simID[isim]),'alpha'] = mae(alpha_obs,alpha[isim],'False')
Metrics.loc[(simID[isim]),'veer'] = mae(veer_obs,veer[isim],'False')
tablecaption = ("**Table 2a. $MAE_{obs}$ with respect to observations**")
display(Markdown(tablecaption))
display(Metrics)
In [36]:
# Errors with respect to reference simulation
refsim = 0
Metrics2 = pd.DataFrame(np.zeros((Nsim, 5)),
columns = ['REWS','Shub','WDhub','alpha','veer'],
index = simID)
Metrics2.loc['units',:] = np.array(['m s-1','m s-1','deg','-','-'])
for isim in range(0,Nsim):
Metrics2.loc[(simID[isim]),'REWS'] = mae(REWS[refsim],REWS[isim],'False')
Metrics2.loc[(simID[isim]),'Shub'] = mae(Shub[refsim],Shub[isim],'False')
Metrics2.loc[(simID[isim]),'WDhub'] = mae(WDhub[refsim],WDhub[isim],'False')
Metrics2.loc[(simID[isim]),'alpha'] = mae(alpha[refsim],alpha[isim],'False')
Metrics2.loc[(simID[isim]),'veer'] = mae(veer[refsim],veer[isim],'False')
tablecaption = ("**Table 2b: $MAE_{ref}$ with respect to **" +simID[refsim])
display(Markdown(tablecaption))
display(Metrics2)
In [35]:
# Metrics vs stability conditions
zLbins = [-20, -0.2, 0.2, 20] # bins
zLbins_label = ['u','n','s'] # stability classes
iz = 1 # 10 m (fluxes at 3m are NaN)
NzL = len(zLbins_label)
zL = zh_obs[iz]/L_obs[iz]
idf = []
for izL in range(0,NzL): idf.append('REWS_'+zLbins_label[izL])
for izL in range(0,NzL): idf.append('Shub_'+zLbins_label[izL])
for izL in range(0,NzL): idf.append('WDhub_'+zLbins_label[izL])
for izL in range(0,NzL): idf.append('alpha_'+zLbins_label[izL])
for izL in range(0,NzL): idf.append('veer_'+zLbins_label[izL])
Metrics3 = pd.DataFrame(np.zeros((Nsim,15)),
index = simID,
columns = idf)
Metrics3.loc['units',:] = np.concatenate((np.repeat('m s-1',NzL),
np.repeat('m s-1',NzL),np.repeat('deg',NzL),
np.repeat('-',NzL),np.repeat('-',NzL)),axis = 0)
for isim in range(0,Nsim):
for izL in range(0,NzL):
i0 = np.where(np.logical_and(zL >= zLbins[izL], zL < zLbins[izL+1]))[0]
Metrics3.loc[(simID[isim]),'REWS_'+zLbins_label[izL]] = mae(REWS[refsim][i0],REWS[isim][i0],'False')
Metrics3.loc[(simID[isim]),'Shub_'+zLbins_label[izL]] = mae(Shub[refsim][i0],Shub[isim][i0],'False')
Metrics3.loc[(simID[isim]),'WDhub_'+zLbins_label[izL]] = mae(WDhub[refsim][i0],WDhub[isim][i0],'False')
Metrics3.loc[(simID[isim]),'alpha_'+zLbins_label[izL]] = mae(alpha[refsim][i0],alpha[isim][i0],'False')
Metrics3.loc[(simID[isim]),'veer_'+zLbins_label[izL]] = mae(veer[refsim][i0],veer[isim][i0],'False')
tablecaption = ("**Table 3: $MAE_{ref}$ for unstable ($u$), neutral ($n$) and stable ($s$) conditions.**")
display(Markdown(tablecaption))
display(Metrics3)
[As per Sanz Rodrigo et al. (2017)]
All of the microscale models produce similar patterns of the diurnal cycle, demonstrating the effectiveness of the offline coupling methodology. Fig. 5 shows contour plots of the evolution of the vertical profile of mean flow quantities computed by the participating models compared to the observations on the top row. The second row corresponds to the reference mesoscale model that is used to derive input forcings and boundary conditions to drive the microscale simulations. Since no additional information is added at microscale, we can assume that this mesoscale simulation is already a good reference for microscale models to verify a correct implementation of the meso-micro methodology. The differences in the mean flow arise from the different ways each model represents turbulence, which is noticed in the $TKE$ contour plots.
The ensemble mean of the WRF simulations used in the sensitivity analysis is plotted in the third row. By ensemble averaging, we obtain a better match with the observations than by using any single WRF simulation – this result was also seen for the LLJ cases discussed in [34]. The ensemble here provides a better prediction of the LLJ with a distinct velocity maximum around midnight instead of a broader double-peak as in the reference WRF simulation.
The microscale models diverge in their estimates of rotor-based quantities of interest. Time series of rotor-based quantities of interest are shown in Fig. 3. The spread of the models for REWS is around 2 $m s^{-1}$ and around 15º in terms of $WD_{hub}$. The spread in terms of wind shear and wind veer is also large specially during nighttime stable conditions. Vertical profiles of wind speed and direction at midnight (Fig. 7) show how the models capture the characteristics of the LLJ. The phase error in the input data dominates the bias in the simulations; this input error cannot be corrected at the microscale by simply changing the turbulence model.
Table 2 summarize the differences between the simulations and observations in terms of $MAE$ with respect to observations (Table 2a, $MAE_{obs}$) and with respect to the reference WRF simulation (Table 2b, $MAE_{ref}$). These are differences integrated over the whole diurnal cycle; therefore mixing all kinds of surface stability and large-scale conditions into a single quantity. We shall focus on the $MAE_{ref}$ to quantify the impact of choosing a different turbulence model at microscale.
WRF-YSU_LES is the closest to the reference, which is to be expected since they are results from different nests of the same simulation. Still, differences are significant of around 0.5 $m s^{-1}$ of wind speed and 4º of wind direction at hub-height. Microscale models increase the error by 0.2-0.7 $m s^{-1}$ of wind speed and up to 7º of wind direction at hub-height with respect to the WRF-YSU_LES results. With respect to observations, all simulations show similar results with a $MAE$ of 1.1-1.6 $m s^{-1}$ of wind speed and up to 14º of wind direction at hub-height. Considering vertical wind speed shear and wind direction veer, SP-Wind and the WRF ensemble produce the closest results to the reference WRF simulation.
Different sets of $k-\epsilon$ constants have been tested (not shown) leading to changes that are within the spread shown in Figure 6, which is also comparable to that observed in the WRF sensitivity analysis when changing the planetary boundary-layer scheme.
Recently, capabilities of EllipSys3D have been extended to cover wall modeled LES of stratified flows and at the same time, a “striped” down 1D version of the code have been made operational [28]. Figures 5 and 6 and Tables 2 and 3 show a general good agreement between URANS based EllipSys1D and EllipSys3D computations. Some minor differences exist thought; a possible cause of them might be related to the fact that the vertical velocity (W) and the advection terms are implicitly assumed to be zero in the EllipSys1D. Further investigations are necessary to confirm this.
Initial LES computations based on coarse grid resolution and very basic Smagorinsky SGS model show that EllipSys3D LES is capable of reproducing the basic LLJ features observed at the Cabauw site, but Tables 2 and 3 MAE comparisons and Figures 2 and 3 indicate that significantly higher resolution and a more advanced approach to SGS turbulence modeling are needed in order to capture all main details relevant for its application in a wind energy context.
Regarding VENTOS®/M wind speed results, both Figures 5 and 7 show that the LLJ magnitude was reasonably well predicted. The diurnal cycle of the simulation wind speed shows over-predictions around 1.4 m s-1, affecting the REWS and Shub error values in Table 2 which are higher than the reference simulation. Good agreement was obtained for the wind direction, shear and veer regarding their integrated error. Analysis of Figure 6 indicates a generalized under-prediction of α for several nocturnal periods, analogous to a higher turbulent shape factor of the boundary-layer. These mismatches happen also with WRF and, despite the VENTOS®/M limitations of its heat-flux boundary condition, the microscale simulation predicts higher values of α and closer to the observations. The results further show a -2.5 K temperature bias that occurs in both WRF and VENTOS®/M results, as well as in the other microscale simulations, which originates from the ERA-Interim input data.
Finally, Table 3 shows the $MAE_{ref}$ is computed for different stability classes filtering with the observed stability parameter $z/L$, where $L$ is the Obukhov length and $z$ = 10 m: unstable ('$u$': $z/L$ < -0.2), neutral ('$n$': -0.2 < $z/L$ < 0.2) and stable ('$s$': $z/L$ > 0.2). Not all the models behave similarly depending on stability. The WRF ensemble REWS is more sensitive in stable conditions. This is also the case for Ellipsys3D_LES probably due to the coarser resolution of the simulation compared to the other LES simulations that do not show this high sensitivity in stable conditions.
In general, it is difficult to extract more meaningful conclusions from Table 2 and Table 3 due to the limited statistical significance of the samples. The overall assessment would be richer if several diurnal cycles from uncorrelated synoptic conditions would have been tested. The ensemble WRF simulations for the same cycle already show significant improvement on mean flow quantities.
[As per Sanz Rodrigo et al. (2017)]
Results of the GABLS3 diurnal cycle benchmark with an emphasis on rotor-relevant values are presented. The main challenge for microscale models was to produce consistent flow fields with respect to the mesoscale model that was used to derive their input forcings. This consistency has been achieved by both LES and URANS models. The spread of the models is significant but of similar magnitude as that shown by WRF using different boundary-layer parameterizations. The input uncertainty coming from the mesoscale, even in relatively ideal conditions, is large and results in $MAE$ of wind speed at hub-height of the order of 1.1-1.6 m s-1 over the whole cycle and and hourly errors of up to 3 m s-1. This is partly mitigated when using an ensemble average of several simulations which also lead to better results in terms of wind shear and wind veer. By ensuring consistency of the microscale models at introducing input forcings we can proceed with further analysis on how RANS and LES models interpret the structure of turbulence in different stability conditions.
This article was produced with funding from the "MesoWake" Marie Curie International Outgoing Fellowship (FP7-PEOPLE-2013-IOF, European Commission’s grant agreement number 624562). It also counts with contributions supported by the New European Wind Atlas (NEWA, FP7-ENERGY.2013.10.1.2 European Commission's grant agreement number 618122). WRF-LES simulations have been produced using high-performance computing resources from the PRACE-MesoWake project in the MareNostrum system, based in Barcelona, Spain. The SP-Wind simulations were performed with computational resources and services provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government – department EWI. JMT’s mesoscale ensemble simulations were conducted using the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. JKL’s efforts were supported by the National Renewable Energy Laboratory under APUP UGA-0-41026-22. NREL is a national laboratory of the U. S. Department of Energy, Office of Energy Efficiency and Renewable Energy, operated by the Alliance for Sustainable Energy, LLC. The benchmark is organized under the umbrella of the IEA-Wind Task 31 "Wakebench" Phase 2 (http://windbench.net/wakebench2).
[1] Sanz Rodrigo J, Chávez Arroyo R-A, Moriarty P, Churchfield M, Kosović B, Réthoré P-E, Hansen KS, Hahmann A, Mirocha JD and Rife D 2016 Mesoscale to microscale wind farm flow modelling and evaluation. WIREs Energy Environ, doi:10.1002/wene.214
[2] Wyngaard JC 2004 Toward Numerical Modeling in the “Terra Incognita”. J. Atmos. Sci. 61: 1816–1826, doi: 10.1175/1520-0469(2004)0611816:TNMITT2.0.CO;2
[3] Holtslag AAM, et al. (2013) Stable atmospheric boundary layers and diurnal cycles—challenges for weather and climate models. Bull. Am. Meteorol. Soc. 94: 1691–1706, doi:10.1175/bams-d-11-00187.1
[4] Baas, P., Bosveld, F.C., Baltink, K. and Holtslag, A.A.M. 2009 A Climatology of Nocturnal Low-Level Jets at Cabauw, J. Appl. Meteorol. Climatol. 48: 1627-1642, doi:10.1175/2009JAMC1965.1, 2009
[5] Bosveld FC, Baas P, vanMeijgaard E, de Bruijn EIF, Steeneveld GJ and Holtslag AAM 2014 The third GABLS intercomparison case for evaluation studies of boundary-layer models, Part A: case selection and set-up. Boundary-Layer Meteorol. 152: 133-156, doi:10.1007/s10546-014-9917-3
[6] Bosveld FC et al. 2014 The third GABLS Intercomparison case for evaluation studies of Boundary-Layer Models Part B: Results and Process Understanding, Boundary–Layer Meteorol. 152: 157-187, doi: 10.1007/s10546-014-9919-1
[7] Sanz Rodrigo J 2017 Windbench/GABLS3 benchmark. http://windbench.net/gabls-3, last accessed in March 2017
[8] Sanz Rodrigo, Churchfield M and Kosovic B 2016 A wind energy benchmark for ABL modeling of a diurnal cycle with a nocturnal low-level jet: GABLS3 revisited. J. Phys. Conf. Ser. 753: 032024, doi:10.1088/1742-6596/753/3/032024
[9] Sanz Rodrigo J, Churchfield M and Kosović B 2017 A methodology for the design and testing of atmospheric boundary layer models for wind energy applications. Wind Energ. Sci. 2: 1-20, doi:10.5194/wes-2-1-2017
[10] Skamarock WC, Klemp JB, Dudhia J, Gill DO, Barker DM, Duda MG, Huang X-Y, Wang W and Powers JG 2008 A description of the advanced research WRF version 3, Technical Note NCAR/TN-475+STR, NCAR, Boulder, CO, June 2008.
[11] Kleczek MA, Steeveneveld GL and Holtslag AAM 2014 Evaluation of the Weather Research and Forecasting Mesoscale Model for GABLS3: Impact on Boundary-Layer Schemes, Boundary Conditions and Spin-Up. Boundary-Layer Meteorol. 152: 213-243, doi:10.1007/s10546-014-9925-3
[12] Hong S-Y and Noh Y 2006 A New Vertical Diffusion Package with an Explicit Treatment of Entrainment Processes. Mon. Wea. Rev., 134: 2318–2341, doi: 10.1175/MWR3199.1
[13] Dee DP et al. 2011 The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. Quart. J. R. Meteorol. Soc. 137: 553-597, doi: 10.1002/qj.828
[14] Wagner R, Cañadillas B, Clifton A, Feeney S, Nygaard N, Martin CSt, Tüxen E and Wagenaar JW 2014 Rotor equivalent wind speed for power curve measurement – comparative exercise for IEA Wind Annex 32. J. Phys. Conf. Ser. 524: 012108, doi:10.1088/1742-6596/524/1/012108
[15] Choukulkar A, Pichugina Y, Clack CTM, Calhoun R, Banta R, Brewer A and Hardesty M 2015 A new formulation for rotor equivalent wind speed for wind resource assessment and wind power forecasting. Wind Energy 19: 1439–1452, doi:10.1002/we.1929
[16] Janjić ZI 1994 The Step-Mountain Eta Coordinate Model: Further Developments of the Convection, Viscous Sublayer, and Turbulence Closure Schemes. Mon. Wea. Rev. 122: 927–945, doi: 10.1175/1520-0493(1994)1220927:TSMECM2.0.CO;2
[17] Nakanishi M and H Niino 2006 An Improved Mellor–Yamada Level-3 Model: Its Numerical Stability and Application to a Regional Prediction of Advection Fog. Boundary-Layer Meteorol 119: 397-407, doi:10.1007/s10546-005-9030-8.
[18] Sukoriansky S, Galperin B and Pero V 2005 Application of a New Spectral Theory of Stably Stratified Turbulence to the Atmospheric Boundary Layer over Sea Ice. Boundary-Layer Meteorol. 117: 231-257. doi:10.1007/s10546-004-6848-4.
[19] Angevine WM 2005 An Integrated Turbulence Scheme for Boundary Layers with Shallow Cumulus Applied to Pollutant Transport. J. Appl. Meteor. 44: 1436–1452, doi: 10.1175/JAM2284.1.
[20] Veiga Rodrigues C, Palma JMLM and Rodrigues ÁH 2016 Atmospheric Flow over a Mountainous Region by a One-Way Coupled Approach Based on Reynolds-Averaged Turbulence Modelling. Boundary-Layer Meteorol.159: 407–437
[21] van de Wiel BJH, Moene AF, Steeneveld GJ, Hartogensis OK and Holtslag AAM. 2007 Predicting the Collapse of Turbulence in Stably Stratified Boundary Layers. Flow Turbulence Combust., 79:251–274
[22] Sogachev A, Kelly M and Leclerc MY 2012 Consistent Two-Equation Closure Modelling for Atmospheric Research: Buoyancy and Vegetation Implementations. Boundary-Layer Meteorol. 145: 307–327, doi:10.1007/s10546-012-9726-5
[23] Deardorff J W 1980 Stratocumulus-capped mixed layers derived from a three dimensional model. Boundary-Layer Meteorol. 18: 495-527, doi:10.1007/BF00119502
[24] Avila M, Folch A, Houzeaux G, Eguzkitza B, Prieto L and Cabezon D 2013 A Parallel CFD Model for Wind Farms. Procedia Comput. Sci. 18: 2157 – 2166, doi:10.1016/j.procs.2013.05.386
[25] Vázquez M et al. 2015 Alya: Multiphysics engineering simulation toward exascale. J. Comput. Scie. 14: 15–27, doi:10.1016/j.jocs.2015.12.007
[26] Michelsen JA 1994 Block structured Multigrid solution of 2D and 3D elliptic PDE’s. Technical Report AFM 94-06, Technical University of Denmark, Department of Fluid Mechanics, May 1994.
[27] Sørensen NN 1995 General Purpose Flow Solver Applied to Flow over Hills. Risø-R-827-(EN), Risø National Laboratory, Roskilde, Denmark, June 1995
[28] van der Laan MP and Sørensen NN 2017 A 1D version of EllipSys. Technical Report, DTU Wind Energy E-0141, March 2017
[29] Meyers J, Meneveau C 2010 Large eddy simulations of large wind-turbine arrays in the atmospheric boundary layer. AIAA Paper No. 2010-827, doi:10.2514/6.2010-827
[30] Munters W, Meneveau C and Meyers J 2016 Turbulent inflow precursor method with time-varying direction for large-eddy simulations and applications to wind farms. Boundary-Layer Meteorol. 159: 305-328, doi:10.1007/s10546-016-0127-z
[31] Allaerts D and Meyers J 2017 Boundary-layer development and gravity waves in conventionally neutral wind farms. J. Fluid Mech. 814: 95-130, doi:10.1017/jfm.2017.11
[32] Moeng C-H 1984 A large-eddy-simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci. 41: 2052-2062, doi:10.1175/1520-0469(1984)0412052:ALESMF2.0.CO;2
[33] Dyer A J 1974 A review of flux-profile relationships. Boundary-Layer Meteorol. 7: 363-372, doi:10.1007/BF00240838
[34] Vanderwende BJ, Lundquist JK, Rhodes ME, Takle ES and Irvin SL 2015 Observing and Simulating the Summertime Low-Level Jet in Central Iowa. Mon. Weater Rev. 143(6): 2319–2336, doi:10.1175/mwr-d-14-00325.1
In [ ]: