atmosphere
ArticleAssessment of Wind over Complex Terrain Considering theEffects of Topography, Atmospheric Stability and Turbine Wakes
Atsushi Yamaguchi, Alireza Tavana and Takeshi Ishihara *Department of Civil Engineering, School of Engineering, The University of Tokyo, Tokyo 113-8656, Japan;
atsushi@bridge.t.u-tokyo.ac.jp (A.Y.); a.tavana@bridge.t.u-tokyo.ac.jp (A.T.)* Correspondence: ishihara@bridge.t.u-tokyo.ac.jp

Abstract: This study proposes a microscale flow model to estimate mean wind speed, fluctuating windspeed and wind direction over complex terrain considering the effects of topography, atmosphericstability, and turbine wakes. Firstly, the effect of topography is considered using Computational
Fluid Dynamics (CFD). Next, a mesoscale model is presented to account for the effect of atmosphericstability. The effect of turbine wakes on the mean and fluctuating wind speeds are then represented
by an advanced wake model. The model is validated using the measurement data of a wind farm
located in the North of Japan. The measured wind data by Lidar at a reference height are horizontally
extrapolated to a nearby met mast hub height and validated by a cup anemometer. Moreover, a novel
averaging method is proposed to calculate a directional equivalent Monin–Obukhov length scale toaccount for the effect of atmospheric stability. Finally, the measured wind data at the reference heightare vertically extrapolated and validated at the lidar location. The predicted mean and fluctuating
wind speeds show good agreement with the measurements.
Keywords: microscale flow model; CFD; mesoscale model; complex terrain; topography; atmosphericstability; turbine wake

1. Introduction
Citation: Yamaguchi, A.; Tavana, A.;
Ishihara, T. Assessment of Wind overComplex Terrain Considering theEffects of Topography, Atmospheric
Stability and Turbine Wakes.
Atmosphere 2024, 15, 723. https://doi.org/10.3390/atmos15060723
Academic Editor: Yubin Li
Received: 27 March 2024
Revised: 10 June 2024
Accepted: 13 June 2024
Published: 17 June 2024

Copyright: © 2024 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms andconditions of the Creative Commons
Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Wind energy is rapidly growing as a renewable energy resource worldwide. In planning prospective wind farms, an accurate prediction of the wind power is crucial to ensure
the economic feasibility of the project. Mesoscale and microscale models have been used
for the assessment of wind. Li et al. [1] used mesoscale model and Computational Fluid
Dynamics (CFD)-based flow model for the local wind prediction. However, accurate windprediction using measurement data is needed for the siting of wind turbine. Tomaszewski
et al. [2] considered the wake of wind farm in the mesoscale model, but the horizontal
resolution was insufficient for the siting of wind turbines. When designing a wind farm,typically, wind velocity measurement is performed at least at one site and they are extrapolated vertically to the hub height of the wind turbine, and horizontally to the planned windturbine positions. According to the WindFarmer theory manual [3], three main sourcesof uncertainty are involved in this process. One source is measurement uncertainty due
to the limitations of the in situ devices; another source is due to the variability of windcondition over time and the third is wind farm modeling. The measurement devices uncertainty can be handled by the procedures explained in IEC61400-12-1 [4]. The variability ofwind condition can also be considered using long-term onsite measurement data [5]. The
most important issue in wind farm modeling is having an accurate wind flow model. On
complex terrain, these vertical and horizontal extrapolations are performed by using flow
models, typically based on CFD. Several research have been performed [6–10] on the flow
models for this horizontal and vertical extrapolation. Bleeg et al. [8] and Meissner et al. [9]used temperature equations in the flow model to consider the effect of atmospheric stability,
but they did not discuss the temperature field. Recently, Liu and Stevens [10] discussed thetemperature field. However, there are still a few problems to be solved.

Atmosphere 2024, 15, 723. https://doi.org/10.3390/atmos15060723

https://www.mdpi.com/journal/atmosphereAtmosphere 2024, 15, 723

2 of 22

The first problem is the cross validation of the flow model for the repowering of theexisting wind farms and development of new wind farms near the existing ones. In thissituation, almost all wind direction is somehow affected by the wake of the existing windturbines and this wake effect must be included in the flow model. In addition, in the sitewhere the atmospheric stability is important, this effect also must be included. When thesurface of the ground is covered by forests, the consideration of forest effect is also needed.
Qian and Ishihara [11] proposed a semi-analytical wake model and carried out a crossvalidation of the wind speed measurement in an existing wind farm considering the wake
effect. However, the atmospheric stability is not considered. Corbett et al. [12] and Uchida
et al. [13] performed the flow simulation by considering the atmospheric stability and forest
effect, but wind turbine wake is not considered.
The second problem is the modeling of atmospheric stability for the flow model.
Strictly speaking, the stability parameter varies with time, and when the downscaling is
performed by using a flow model, all the combinations of wind direction and stability
must be considered, which is computationally expensive. To reduce the computationalcost, Durán et al. [14] proposed to use the representative stability for each wind direction
by taking the harmonic average of Monin–Obukhov Length scale. However, the effect ofMonin–Obukhov length scale on the wind speed is highly nonlinear, and it is not clear
whether harmonic averaging is applicable or not.
The last problem is that the effects of topography, atmospheric stability, and wakemodel on the vertical profile of mean wind speed and turbulence intensity need to beconsidered at the same time and a framework to account for these effects is necessary.
In this study, a framework for a microscale flow model to account for the effects of
topography, atmospheric stability, and wake from other turbines is proposed. In Section 2,
firstly, a CFD model to evaluate the effect of topography is presented. A stability correction
factor is then proposed to consider the effect of atmospheric stability using reverse Monin–
Obukhov length (1/L) data extracted from a mesoscale model. Finally, an advanced wakemodel is introduced to account for the effects of turbine wakes. The predicted mean andfluctuating wind speeds are validated with the measurements by a vertical Lidar and a Met
mast in Section 3. The conclusions are summarized in Section 4.
2. Methodology
The overview of the proposed approach is shown in Figure 1. Based on the onsite
measurement data, horizontal and vertical extrapolation of wind speed is performed byconsidering the effects of topography, atmospheric stability, and wind turbine wakes as
described in Sections 2.1–2.3.
First, CFD model is used to calculate speed up ratio and wind direction change causedby the topography and surface roughness. The details of the CFD model are described in
obs , fluctuating wind speed, σobs , and wind
Section 2.1. If measured mean wind speed, Urefrefobs are available at one point, by using terrain correction factors, these variables
direction, θrefat any point ( x, y, z) can be calculated.
terrainobs
terrainobs
Utargetx, θref(1)(x, t) = Uref(t) × CUterrainobs
obs
σtarget(x, t) = σref(t) × Cσterrain x, θref
terrainobs
obsθtarget(x, t) = θref(t) + ∆θ terrain x, θref

(2)
(3)

terrain , Cterrain are the terrain correction factors for the mean and fluctuatingwhere CUσwind speeds, respectively, and ∆θ terrain is the correction term for the wind direction andx = ( x, y, z) is the position of the site where the wind is calculated.

Atmosphere 2024, 15, x FOR PEER REVIEW

3 of 24

where 𝐶𝑈t𝑒𝑟𝑟𝑎𝑖𝑛 , 𝐶𝜎t𝑒𝑟𝑟𝑎𝑖𝑛 are the terrain correction factors for the mean and fluctuating
3 of 22wind speeds, respectively, and ∆𝜃 t𝑒𝑟𝑟𝑎𝑖𝑛 is the correction term for the wind direction and𝐱 = (𝑥, 𝑦, 𝑧) is the position of the site where the wind is calculated.

Atmosphere 2024, 15, 723

Mesoscale model

CFD simulation

Onsite measurement

obs
Mean wind speed 𝑈ref(𝑡)

Terrain correction factor for mean wind speed 𝐶

obs
Fluctuating wind speed 𝜎ref(𝑡)
obs
Wind direction 𝜃ref (𝑡)Terrain correction factor for fluctuating wind speed 𝐶
Terrain correction term for wind direction 𝜃(𝐱, 𝜃)

Monin-Obukhov length scale 𝐿(𝑥, 𝑦, 𝑡)(𝐱, 𝜃)(𝐱, 𝜃)

terrainTerrain corrected mean wind speed 𝑈target(𝐱, 𝑡)

stabilty

Stability correction factor𝐶𝑈(𝐱, 𝜃)

terrainTerrain corrected fluctuating wind speed 𝜎target(𝐱, 𝑡)
terrainTerrain corrected Wind direction 𝜃target(𝐱, 𝑡)

initialWake free mean wind speed 𝑈target(𝐱, 𝑡)
initialWake free fluctuating wind speed 𝜎target(𝐱, 𝑡)
initialWake free wind direction 𝜃target(𝐱, 𝑡)

Analytical wake model

pred

Predicted wind speed 𝑈target (𝐱, 𝑡)
pred

Predicted fluctuating wind speed 𝜎target (𝐱, 𝑡)
pred

Predicted wind direction 𝜃target (𝐱, 𝑡)

Figure 1. A framework of the microscale flow model.
Figure 1. A framework of the microscale flow model.

Then, thethe Monin–Obukhov
extracted fromfrom thethe mesoscale
mesoscale simulation
simulation
Then,
Monin–Obukhov length
length scale
scale L𝐿 isisextracted
stabilityobs
stability
and
used
to
calculatethestability
correction
factor
𝐶
(𝐱,
𝜃
)
which
is
also
a funcobsref which is also a function
and used to calculate the stability correction factor
CU 𝑈x, θrefobs
tion of the point (𝐱) and wind direction𝜃ref . The detail of this correction factor is deobs . Theof the point (x) and wind direction θrefdetail of this correction factor is described in
scribedin
Section
2.2.
AccordingtoStiperski
andtheCalaf[15],
effectonof fluctuatingstability onwind
flucSection 2.2. According to Stiperski and Calaf [15],effectofthe
stability
tuating
wind
speed
is
small
in
neutral,stable,andmoderately
unstablecases.
So,in
this
speed is small in neutral, stable, and moderately unstable cases. So, in this study, the effectstudy,
theeffectofstabilityonthefluctuating
wind
speed
is
neglected.
of stability on the fluctuating wind speed is neglected.
stabilityobs 
initialterrain𝑈target(𝐱, 𝑡) = 𝑈target(𝐱, 𝑡) × 𝐶𝑈 stability(𝐱,𝜃ref)initialterrainobs
Utarget (x, t) = Utarget (x, t) × CU
x, θref(4)(4)

initial (𝐱,terrain (𝐱,𝜎target𝑡) = 𝜎target𝑡)(5)
initialterrainσtarget (x, t) = σtarget(5)(x, t)
initial (𝐱,terrain𝜃target
𝑡)
=
𝜃
(𝐱,
𝑡)
(6)
initialterrainθtargetθtarget(6)(x, t) =target(x, t)
Finally, an advanced wake model is applied to calculate the velocity deficit
Finally, an advanced wake model is applied to calculate the velocity deficit ∆U ( x, y, z, t),∆𝑈(𝑥, 𝑦, 𝑧, 𝑡), and added turbulence ∆𝜎(𝑥, 𝑦, 𝑧, 𝑡) from wind turbines. The detail of theand added turbulence ∆σ( x, y, z, t) from wind turbines. The detail of the wake model iswake model is described in Section 2.3.
described in Section 2.3.
predinitial (𝐱,𝑈target (𝐱, 𝑡) = initial
𝑈target𝑡) − ∆𝑈(𝐱, 𝑡)(7)
pred
Utarget (x, t) = Utarget (x, t) − ∆U (x, t)(7)

q initial 22 (𝐱, 𝑡)2 (𝐱, 𝑡) + ∆𝜎√𝜎initialtarget(x,σtargett) + ∆σ2 (x,
t)
pred
pred
𝐼targetItargett) 𝑡)==(x,(𝐱,
pred
pred (𝐱, 𝑡)
U𝑈(x, t)
targettargetpredinitialθtarget (x, t) = θtarget(x, t)

(8)
(8)
(9)

Atmosphere 2024, 15, 723

4 of 22

pred

pred

pred

where Utarget (x, t), Itarget (x, t) and θtarget (x, t) are the predicted mean wind speed, turbulence intensity and wind direction, respectively. The effect of multiple turbines issuperposed to calculate the final values as proposed by Qian and Ishihara [16].
By substituting Equations (7) and (4) into Equation (1), the following equation for themean wind speed can be obtained.




pred
stability
obs
terrainobs
obs
Utarget (x, t) = Urefx, θref× CU
x, θref− ∆U (x, t)(10)(t) × CU
This equation is based on dimensional analysis and the effects of terrain and atmoterrain
spheric stability are accounted for by multiplying the non-dimensional coefficients CU
stability

and CU. This paper aims to propose a method to identify these coefficients.

2.1. Computational Fluid Dynamics
To consider the effect of local topography and surface roughness, CFD simulations
can be used. The momentum is mainly driven by large eddies. In this study, the Large
Eddy Simulation (LES) model is used in which large eddies are resolved and the eddies
smaller than the grid spacing are parameterized. To calculate these sub-grid scale (SGS)
stresses, the Smagorinsky–Lilly model is used. The conservation of mass and momentum
are shown in Equations (11) and (12). ANSYS Fluent is used as the CFD solver.
∂ρuej
=0
∂x j
∂ρuei uej
∂ρuei
∂e
p
∂
+
=−
+
∂t
∂x j
∂xi
∂x j

(11)
∂ue
µ i
∂x j

!

+

∂τij
+ f u,i∂x j(12)where uei is the filtered wind speed in the xi direction. ρ is the air density, pe is the filteredpressure and µ is the molecular viscosity of the air.
τij is the SGS stresses which is modeled as follows:1τij = 2µt Seij + τkk δij
3
!
∂uej
e
1
∂
u
i
Seij =+
2 ∂x j
∂xi

(13)

(14)

where µt denotes the SGS turbulent viscosity, Seij is the rate of strain tensor for the resolved
scale and δij is the Kronecker delta.
Smagorinsky–Lilly is used for the SGS turbulent velocity:µt = ρL2s Se = ρLsq

2Seij Seij
Ls = min κd, Cs V 1/3

(15)
(16)

where Ls denotes the mixing length for the sub-grid scales, κ is the von Karman constant,i.e., 0.42, d is the distance to the closest wall and V is the volume of the computational cell.
In this study Cs , the Smagorinsky constant is set to be 0.032.
f u,i is the source term added to the momentum equation to account for the effect of
obstacles using the canopy model as shown in Ishihara et al. [17]. Fluid force on bluff body,
f u,i , is derived in Equation (17).
f u,i = −

Fu,i
1
1 q
CD
= ρC f γ0 u2j ui , C f =
Vgrid
2
l0
(1 − γ0 )2(17)

Atmosphere 2024, 15, 723

5 of 22

where c f is a kind of drag coefficient representative, l0 is the representative length scaleof the obstacle and γ0 is the packing density. The values of c f , l0 and γ0 need to be
set according to different types of obstacles. To consider the effect of vegetation, zerodisplacement plane approach as shown by Oke [18] is adopted, i.e., the vertical axis isshifted by the zero-plane displacement height.
Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm is used
for the pressure velocity coupling, as shown by Ferziger and Peric [19]. Finite volume
method and unstructured collocated mesh are used for the present simulations. The second
order central difference scheme is used for the convective and viscosity term, and thesecond order implicit scheme for the unsteady term. The CFD simulation is executed in 16
representative direction sectors.
The CFD results can be used to obtain the mean and fluctuating components of windspeed, and wind direction at any point inside the domain using reference point data bycorrecting the effect of terrain and roughness in each direction sector. The correction factor
for each of these variables is defined in Equations (18)–(20).
eterrain (x, nθ ) =CUterrain (x, n )Utargetθ
terrain (xUrefref, , nθ )eσterrain (x, nθ ) =Cterrain (x, n )σtargetθ
terrain (xσrefref, nθ )terrainterrainf terrain (x, nθ ) = θtarget∆θ(xref , nθ )(x, nθ ) − θref

(18)

(19)
(20)

in which σterrain is the fluctuating wind speed, U terrain is the mean wind speed and θ terrain
is the wind direction obtained from the CFD results and nθ is the index number of diobs , are obtained using a linear
rection sector. The correction factors for any direction, θrefinterpolation.
obs is obtained as shown in Equation (21).
First, the direction sector of θref
nθ =



obs +11.25
 θref22.5

obs < 348.75
θref

0

(21)

else

The linear interpolation is then performed as follows.

terrainobs
terrainterrainterrainCUx, θref= CeU(x, nθ ) + a CeU(x, nθ + 1) − CeU(x, nθ )




obs
Cσterrain x, θref= Ceσterrain (x, nθ ) + a Ceσterrain (x, nθ + 1) − Ceσterrain (x, nθ ) terrain
f terrain (x, nθ ) + a ∆θff terrain (x, nθ )∆θ terrain x, θ obs = ∆θ(x, nθ + 1) − ∆θ
ref

(22)
(23)
(24)

where the proportional coefficient a is defined in Equation (24).
a=obs − θ terrain ( n )
θref
θ
ref
terrain ( n + 1) − θ terrain ( n )
θref
θ
θ
ref

(25)

2.2. Mesoscale Model and Stability Correction
As a mesoscale model, the Weather Research and Forecasting (WRF, 3.4.1) modelingsystem, a state-of-the-art software for atmospheric research and forecasting applications,is used. WRF has proven to be a reliable choice in a multitude of scales from hundredsof meters up to thousands of kilometers [20]. The governing equations in the Advanced
Research WRF (ARW) solver are briefly explained below. The parameters used in theequations are listed in Table 1.

Atmosphere 2024, 15, 723

6 of 22

Table 1. List of parameters in WRF.
Types

Parameter

Definition

Prognostic variablesv(u, v, w)
ϕ = gzθmq v , q c , qr , qi

Covariant velocities
Geopotential
Moist potential temperature
mixing ratios for water vapor, cloud, rain, ice

Diagnostic variables

p
Pdh
Pht
αd = 1/ρθµdPressure
Hydrostatic component of pressure for dry air
Ph along the top boundary for dry air
Reciprocal of dry atmospheric density
Potential temperaturemass of dry air in the column

Constantsgγ = C p /Cv
p0
Rd
Rv

Gravitational acceleration
Specific heat ratio
Reference pressure 105 Pa
Gas constant for dry air
Gas constant for water vapor

Predefined variables

η

Vertical coordinate

Forcing terms

FU , FV FW , FΘ , FQm

External forces

The ARW equations are formulated using a terrain-following hydrostatic pressure
vertical coordinate denoted by η.
η=Pdh − Pdht
µd

(26)

where µd is the mass of dry air in the column and Pdh and Pdht are hydrostatic pressure of thedry atmosphere and the hydrostatic pressure at the top of the dry atmosphere, respectively.
The coupled variables are defined as V = µd v and Θ = µd θ. The flux equations are: ∂p
α ∂p ∂ϕ∂U+ (∇·Vu) + µd α +
= FU
(27)
∂t
∂x
αd ∂η ∂x ∂V∂p
α ∂p ∂ϕ+ (∇·Vv) + µd α +
= FV
(28)
∂t
∂y
αd ∂η ∂y
 

∂W
α ∂p
+ (∇·Vw) + g− µd = FW
(29)
∂t
αd ∂η
The governing thermodynamics equation, conservation of mass, geopotential time
derivative and the equation of state are, respectively, shown in Equations (30)–(33).
∂Θ+∇·Vθ = FΘ
∂t

(30)

∂µd
+ ∇·V = 0∂t∂Φ1
+ [(V·∇ϕ) − gW ] = 0
∂t
µd


Rd θm γp = p0p0 α d(31)∂Qm+ (∇·Vqm ) = FQm
∂t

(32)
(33)
(34)

Atmosphere 2024, 15, 723

7 of 22

where α is the inverse density taking the whole parcel
α = αd (1 + qv + qc + qr +qi )−1andQm = µd qm ; qm = qv , qc , qi .

intoaccountv
Additionally, θm = θ 1 + R
Rd qv ≈ θ (1 + 1.61qv ).
In the mesoscale model, for assimilating the reanalysis data, nudging with reanalysis
data is used to keep the simulations close to reanalysis over the course of simulation asshown in Equation (35) [21].

∂θ= F (θ ) + Gθ θb0 − θ
∂t

(35)

where θ is the prognostic variable, F (θ ) represents the normal tendency terms due tophysics, the nudging parameter, Gθ is the timescale controlling the nudging strength and θb0
is the time and space interpolated analysis field value towards which the nudging relaxes
the solution.
In most CFD simulations, the assumption is that the atmosphere has a neutral atmospheric stability. However, this might not be true in many cases. Therefore, the atmosphericstability effect needs to be considered in case such effect exists [9]. This effect can be
included using mesoscale simulation results. The wind speed profile can be corrected using
the following formula [22–26].
  
 z u∗zu(z) =
ln
−ψ
(36)
κ
z0
L
where u∗ is the friction velocity, κ is the von Karman constant, z0 is the surface roughness
and L is the Monin–Obukhov
 length scale, which is extracted from the mesoscale simulation
results; the function ψ Lz is defined in Equation (37).

ψ

z
L

=


−15.14
 −0.35 z


z
z
5

L − 0.75×5
−
0.7
−
0.75
−
e

L
L
0.35
0.35




−5 z
L

0
2








2 ln 1+2 x + ln 1+2x − 2 tan−1 ( x ) + π2





1/4

x = 1 − 16 Lz


7 ≤ Lz
0.5 < Lz < 7
0 < Lz ≤ 0.5
z
L =0

z
L <0

(37)

The CFD results can be corrected to include the effect of stability at time t. This value
is defined as the height correction ratio of stratified profile to neutral profile over flat
terrain as:
i
h  

stability

CU_t(x, t) =z
ln zz −ψ L( x,y,th z 0  z ) irefln zref −ψ L( x,y,t)0h  iln zz
h  z 0 iln zref(38)

0

The idea for Equation (37) is that the effect of topography and the effect of atmostability
spheric stability can be separated, and the stability correction factor CU_t
(x, t) is definedstability

assuming the terrain is flat. The value CU_t
(x, t) becomes unity for the same height(z = zref ) or neutral (ψ = 0). This approach requires time history analysis and it is difficult
to discuss the directional vertical profile. As discussed earlier, Duran et al. [14] proposed anequivalent value of stability correction factor by taking the harmonic average of L by using
the data where wind speed is larger than 3 m/s. In this study, an alternative approachAtmosphere 2024, 15, 723

8 of 22is proposed to calculate the directional stability correction factor. The proposed factor is
calculated using Equation (39).

estability (x, nθ ) =
C
U

h  
i
ln zz −ψeq (x,nθ )h  z 0iln zref −ψeq ( x,y,zref ,nθ )0h  iln zz
h  z 0 iln zref(39)

0where ψeq (x, nθ ) can be calculated as shown in Equation (39).

ψeq (x, nθ ) =obs ( t ) ψ∑ Urefz
L( x,y,t)(40)

obs ( t )∑ Urefobs ( t ) is the observed wind speed in that direction sector at time t. Once ψ
Uref
eq is
calculated for each direction sector, the equivalent value of 1/Leq for this wind direction
can be calculated inversely which can explain the effect of stability for this wind directionand expressed as:ψ−1 ψeq (x, nθ )1(41)=
Leq ( x, y, nθ )zstability
obs at any direction θ obs , the C
estability (x, nθ ) is linearly interpoTo obtain CUx, θrefU
reflated in the similar way as explained in Section 2.1.

2.3. Turbine Wake ModelIn the downstream of turbines, the wind is affected by turbine wakes. In these cases,
the effect of the wake needs to be considered and corrected. Turbine wake is a topic thathas been widely investigated and analytical models have been developed to account for it.
In this study, the effect of wake on the mean and fluctuating components of wind speed
is calculated using the Ishihara–Qian wake model; the summary is shown in Table 2 [16].
This analytical Gaussian-based model provides the three-dimensional wake characteristics
including wake width, velocity deficit and added turbulence.
Table 2. Summary of the Ishihara–Qian wake model.
Wake Model

Formulaσ∗ x∗
D = k D + ε , Dw = 4

Wake width∆U ( x,y,z)= nUh

Velocity Deficitr=

Added Turbulence

Parameters1−2d+e· Dx + f (1+ Dx )

(
∆Ia =



1
a+b· Dx +c(1+ Dx )

q

−2

√

1.07

k∗ = 0.11C ′ T Ia0.20
0.17
T − 0.25Ia
C ′ T = CT (Uh , γ)cos γ, Ia > 0.03ε∗ = 0.23C ′2 ln 2σo2 exp



2

r
− 2σ
2

−0.75

a = 0.93CC ′ T
Ia0.17
0.6 0.2′
b = 0.42C T Ia
−0.25 −0.7
c = 0.15C ′ T
Ia
−1.2
′
d = 2.3C T , e = 1.0Ia0.1

( y − y d )2 − ( z − H )2∆σ( x,y,z)
=
Uh




k1 exp −(r − D/2)
2σ2

2



r + D/2+ k2 exp − ( 2σ2 )

2



else

 0
y
Ia sin2 π HD−z cos2 π D (0 ≤ z < H, |y| ≤ D )− ∆Ia

f = 0.7CT′−3.2 Ia−0.45
 rcos2 π2 Dr − 0.5
D ≤ 0.5
k1
r
> 0.5
1
 2 π r D
rcos 2 D + 0.5
≤
0.5
D
k1
r
0
D > 0.5


where D is the rotor diameter, CT is the thrust coefficient, Uh is the rotor onset wind speed and Ia is the rotor onset
turbulence intensity.

As the turbines experience a non-uniform flow caused by the complex terrain or thewake of wind turbines further upstream, it is first required to calculate the rotor onset
wind speed Uh and turbulence Ia . The rotor onset wind speed Uh is calculated by directlyAtmosphere 2024, 15, 723

9 of 22

performing a geometric averaging of wind speed U over the rotor, while rotor onset Ia iscalculated by the root mean of squares of streamwise turbulence standard deviation σu
over the rotor divided by the rotor onset wind speed Uh [16].
Uh =

1
A

Z

1
Ia =
AUh

rotor

rZ

U ( x, y, z)dA(42)

σu2 ( x, y, z)dA

(43)

Using these values, the effect of wake on the flow field, ∆U and ∆σ, can be calculated
for each wind turbine.
3. Results and Discussion
The on-site measurement data from a wind farm located in the north of Japan are
described in Section 3.1 and the proposed microscale flow model are validated in Section 3.2.
The modeling of atmospheric stability and its effect on vertical extrapolation is discussed
in Section 3.3.
3.1. Description of the Site and Computational Setting
To validate the proposed methodology, the data from Iwaya wind farm located in
the northern Japan is used in this study. The location of this wind farm can be found in
Figure 2a. A Met mast and a vertical profiler lidar, WindCube V2, are installed in the windfarm as shown in Figure 2b. At the mast, wind speed and wind direction are measured at
the hub-height, which is 68 m above ground level and the lidar measures wind velocity
every 10 m between 40 m and 230 m above ground level. Seven months of measurement
data from 23 November 2018 to 23 June 2019 are used in this study. The wind speed
measurement by using vertical lidar on complex topography includes errors due to theassumption of homogenous flow in their wind vector reconstruction as shown in [27]. Thiserror is estimated and compensated using the CFD simulation as explained in Appendix A.
Atmosphere 2024, 15, x FOR PEER REVIEW
of 24
The position of the wind turbines and the location of the lidar and Met mast sites are 10
shown
in Figure 2b.

Y
Metmast
Lidar

x

100 m

(a)

(b)

Figure 2.
2. Iwaya
Iwaya wind
wind farm.
farm. (a)(a) Location;
Location; (b)(b) Turbines,
Turbines,met
metmastmastandandlidarlidarlayout.
layout.
FigureNNW

At the Met mast and lidar locations, several wind direction sectors are affected byNturbinewakes. Figure 3 shows the alignmentNNWof turbines
relativeN
NNE to the lidar and met mast.
NNE

NW

NE

NW

ENE

WNW
Metmast

WNW

NE

ENE

Atmosphere 2024, 15, 723

10 of 22(a)(b) exist within 10D, where D = 62 m
These wake-affected directions in which wind turbines
is therotor diameterin thiswind farm,shaded met
in thefigure.
Figure2. Iwayawind farm.
(a) Location;(b)areTurbines,
mast
and lidar layout.

Atmosphere 2024, 15, x FOR PEER REVIEW

11 of 24mesh size of which is approximately 2000 m by 2000 m, is covered by using a constanthorizontal resolution of 30 m, and the mesh size near the center (Figure 5b) is 7.5 m, basedon Qian and Ishihara [29]. They investigated the mesh size dependency and concluded
that there is a difference between 16 m and 8 m, but the difference between 8 m and 4 m
is negligible. In this study, the turbulence intensity is calculated based on the unsteadycomponent only since the contribution from the subgrid turbulence is negligible for meshsizes less than 8 m. Outside this area, the horizontal mesh size is stretched with a constantratio of 1.15 up to 240 m. The computational target domain size of 8000 m by 8000 m (Fig(a) is defined so that almost all the important topographic(b)
ure 5c)characteristics are included.
Furthermore,
additional
domainto the upstream
of(a)
the
computationaltargetWake-affecteddomain and
Figure3. Theplacement
ofoftherelativetoto
Met
mast
and
Lidar.
Figure3. Theplacement
thewindwindturbines
turbines relative(a) Metmastand(b)(b)
Lidar.
Wake-affectedbufferzones
for the boundary conditions are added making the total computational dodirections
areareshaded.
directions
shaded.
main size of 30,000 m in wind direction and 20,000 m in cross wind direction (Figure 5a).
The minimum
verticalmeshsize isdistribution
5 m at the groundandspeedstretches
constant
ratio
windroseandandfrequency
distribution
ofofthewind
at lidar160them are
TheThewindrosefrequency
the
wind
speed
atwithlidar
160
mshownare shownin
Figure
4.
of
1.1
and
the
top
boundaryislocated
10,000m
above
sea
level.
in Figure 4.
The computation domain of the CFD simulation is shown in Figure 5a and the meshis shown in Figure 5b. The minimum mesh domain size of 2000 m is based on 30 times of
the hub height (68 m). Rotating mesh is used for each wind direction and horizontal dimensions include the upstream additional domain as described in Ishihara et al. [28]. A
typical example of the computational mesh for the westerly wind with the topography
and surface roughness map is shown in Figure 5c. Terrain following coordinate is used
and the boundary conditions are set as Qian and Ishihara [29]. As the typical height of thetrees in the area is 10 m, a zero-displacement plane of 7 m is considered on the bottom
surface. To generate the input data, two rows of blocks are placed in the upstream as explained in Qian and Ishihara [16]. The timestep of the simulation is 1 s. The mesh settingof the CFD simulation is shown in Table 3. The wind farm area shown in Figure 2b, themesh size of which is approximately 2000 m by 2000 m, is covered by using a constanthorizontal resolution of 30 m, and the mesh size near the center (Figure 5b) is 7.5 m, based

(a)

(b)

Figure4. Windwindfarm
farm
lidar.
(a) Windrose;(b) Frequency
distribution
the windFigure4. WindatatIwaya
Iwaya windlidar.
(a) Windrose; (b)Frequency
distributionof the windof
speed.
speed.
The computation domain of the CFD simulation is shown in Figure 5a and the meshis shown in Figure 5b. The minimum mesh domain size of 2000 m is based on 30 timesof the hub height (68 m). Rotating mesh is used for each wind direction and horizontal
dimensions include the upstream additional domain as described in Ishihara et al. [28]. A
typical example of the computational mesh for the westerly wind with the topography andMet mast
surface roughness map is shown in Figure 5c. Terrain following coordinate
is used and theLidarratio of 1.15 up to 240 m. The computational target domain size of 8000 m by 8000 m (Figure 5c) is defined so that almost all the important topographic characteristics are included.
Furthermore, additional domain to the upstream of the computational target domain andbuffer zones for the boundary conditions are added making the total computational do11 of
22main size of 30,000 m in wind direction and 20,000 m in cross wind direction (Figure
5a).
The minimum vertical mesh size is 5 m at the ground and stretches with the constant ratioof 1.1 and the top boundary is located 10,000 m above sea level.
boundary conditions are set as Qian and Ishihara [29]. As the typical height of the trees
in the area is 10 m, a zero-displacement plane of 7 m is considered on the bottom surface.
To generate the input data, two rows of blocks are placed in the upstream as explainedin Qian and Ishihara [16]. The timestep of the simulation is 1 s. The mesh setting of the
CFD simulation is shown in Table 3. The wind farm area shown in Figure 2b, the mesh sizeof which is approximately 2000 m by 2000 m, is covered by using a constant horizontal
resolution of 30 m, and the mesh size near the center (Figure 5b) is 7.5 m, based on Qianand Ishihara [29]. They investigated the mesh size dependency and concluded that there isa difference between 16 m and 8 m, but the difference between 8 m and 4 m is negligible.
In this study, the turbulence intensity is calculated based on the unsteady component only
since the contribution from the subgrid turbulence is negligible for mesh sizes less than
8 m. Outside this area, the horizontal mesh size is stretched with a constant ratio of 1.15 up
to 240 m. The computational target domain size of 8000 m by 8000 m (Figure 5c) is definedso that almost all the important topographic characteristics are included. Furthermore,
additional domain to the upstream of the computational target domain and buffer zones
for the boundary conditions are added making the total computational domain size of
30,000(Figure 5a). The minimum(a) m in wind direction and 20,000 m in cross wind direction(b)vertical mesh size is 5 m at the ground and stretches with the constant ratio of 1.1 and the
Figure 4. Wind at Iwaya wind farm lidar. (a) Wind rose; (b) Frequency distribution of the wind
top boundary is located 10,000 m above sea level.
speed.

Atmosphere 2024, 15, 723

(a)

(b)

(c)

Figure5. CFD
Figure 5.
CFD simulation.
simulation. (a)(a) Domain
Domain andand boundary
boundary conditions;conditions; (b)(b) Horizontal
Horizontal meshmesh onon thethe bottom
bottomsurface; (c) Typical mesh and the location of the lidar and the met mast for westerly wind direction.
surface; (c) Typical mesh and the location of the lidar and the met mast for westerly wind direction.
Table 3. CFD simulation domain and mesh setting.
Table 3. CFD simulation domain and mesh setting.

Domain and Mesh Setting
Domain and Mesh Setting
Minimum horizontal mesh size
7.5 m
Minimumhorizontal
mesh sizemesh size
7.5 m 240 m
Maximumhorizontal
Maximum horizontal mesh size
240 m
Minimum mesh domain size
800 m (Diameter)
Minimum mesh domain size
800 m (Diameter)
Horizontal
mesh stretching
ratio
1.15
Horizontal
mesh stretching
ratio
1.15
Minimum
vertical verticalmesh size
5m
Minimummesh size
5m
Verticalmesh stretching
ratio
1.1
Verticalmesh stretching
ratio
1.1
For the mesoscale simulation, WRF Version 3.4 is used in periods of one month from
November 2018 to June 2019 when the measurement data is available. The simulations are
conducted with a spin-up period of two days for each month. The computational domain
and the model configuration used in the WRF simulations are shown in Figure 6 and Table 4.
Three nested domains are set around Iwaya wind farm in Aomori. The outer domain has
an 18-km resolution, the inner domain has a 6-km resolution, and the innermost domain
has a resolution of 2 km. Forty-five layers are employed in the vertical direction. Terrain
data provided by the Geographical survey institute of Japan and land use data provided
by the ministry of land, infrastructure, transport, and tourism are used. The Operational
Sea Surface Temperature, and Ice Analysis (OSTIA) dataset is used [30] as the sea surface
temperature.

Atmosphere 2024, 15, 723

has an 18-km resolution, the inner domain has a 6-km resolution, and the innermost domain has a resolution of 2 km. Forty-five layers are employed in the vertical direction.
Terrain data provided by the Geographical survey institute of Japan and land use data
provided by the ministry of land, infrastructure, transport, and tourism are used. The Op22
erational Sea Surface Temperature, and Ice Analysis (OSTIA) dataset is used [30]12asof thesea surface temperature.

Figure6.6.Computation
Computationdomain
domainconfiguration
configurationfor
forWRF
WRFsimulations.
simulations.
FigureTable 4. Description of mesoscale model and computational conditions.
Table 4. Description of mesoscale model and computational conditions.

Domain 1

Simulation period
Simulation
periodCenter point
Center point

Horizontal resolution
Horizontal resolution

Vertical resolution
Verticalresolution
Time
step
Spin-upTime
step
Boundary
condition
Spin-up
Boundary condition

Physics option
Physics option

Domain 2
Domain 3
Domain 2
Domain 3
November 2018–June 2019
November
2018–June
45°N
143°E 2019◦
N 143◦ E
18 km
645km
2 km
(100 × 100
grids)(130 × 1546 km
grids)(100 × 100
18 km
2 km grids)(100 × 100 grids)(130(Surface
× 154 grids)(100 × 100 grids)
45 levelsto 50 hPa)
45 levels24(Surface
72 sec.
sec. to 50 hPa)
8 sec.
Two
72 s
24 dayss
8s
Era 5 0.25°Two× 0.25°dayshourly [31]Ferrier (newEta)
microphysics[32]◦
Era 5 0.25 × 0.25◦ hourly scheme[31]
Rapid radiative transfer model [33]Ferrier (new Eta) microphysics scheme [32]Dudhia scheme[34] [33]Rapid radiativetransfer model
Unified Noah
landsurfaceDudhia scheme [34]model [35]UnifiedNoahmodel2.5[35]
Mellor–Yamada–Janjicland(Eta)surface
TKE level
scheme [36]Mellor–Yamada–Janjic
(Eta)
TKE
level
2.5
scheme[36]
Betts–Miller–Janic scheme (Except for domain 3)[37]
Domain 1

Betts–Miller–Janic scheme (Except for domain 3) [37]In this study, Gθ is set as 3.0 × 10−4 and grid nudging to the reanalysis data (ERA5)
is used. The nudging parameters are the wind velocity and potential temperature. The
nudging scheme proposed by Kikuchi et al. [38] is used in which nudging is applied at all
the vertical levels in domain 1 and only above 6 km in domains 2 and 3.
3.2. Cross Validation of the Flow Model
Cross validation of the two available measurement data is performed by using the
proposed framework. The lidar-measured wind speed (160 m above ground level) andwind direction (70 m above ground level) are used as reference data and wind speed anddirection at the met mast are predicted by using the proposed method. Very low windspeeds cause large uncertainties and are not important in power generation. Therefore,
wind speeds lower than 4 m/s [39] are filtered out.

Atmosphere 2024, 15, 723

13 of 22

The predicted mean wind speed at the met mast at the hub height is plotted as afunction of wind direction and compared with the measurement data in Figure 7a. The redline shows the case where only the effect of terrain is considered. This approach causes
overestimation of the mean wind speed for some wind directions due to turbine wakes(i.e., N, NNE, NE, S, and SW). When the effect of atmospheric stability is considered (blueline), the overestimation is partly improved (e.g., 180 to 225 degrees of wind direction, i.e.,
S to SW). Consideration of the turbine wake improves the predicted mean wind speeds
further (green line), especially in the wind directions when the met mast is in some of
the turbine wakes (e.g., 0, 22.5, 45, 180, 202.5 and 225 degrees of wind direction (i.e., N,
NNE, NE, S, SSW and SW). For the wind direction 202.5 degrees (i.e., SSW direction)
both atmospheric stability and turbine wake are important. The result of the mesoscalemodel WRF simulation downscaled by using the microscale model is also shown in the
purple dotted line. Obviously, the mesoscale model has larger errors, although all of the
microscale effects are included. The frequency distribution of the wind speed at the met
mast at the hub height is also shown in Figure 7b. The overestimation of the wind speed isclearly improved by considering the effects of atmospheric stability and turbine wake. The
prediction accuracy of the wind direction is also improved by using local measurement
Atmosphere 2024, 15, x FOR PEER REVIEW
14 of 24
data compared to WRF with the microscale model as shown in Figure 7c. It shouldbenoted that wake model and stability model does not change the wind direction as shownin Equations (1)–(9), and predicted frequency distribution of wind directions are identical
for the three models.

Mean wind speed [m/s]20
Measurement
TerrainTerrain+Atmospheric stability
Terrain+Atmospheric stability+Wake
WRF+Microscale Model

15

10

5

0
0

90

180

270

360

wind direction [degree]

(b)

(a)

(c)
Figure7. Validationof windspeed
at thematmatmast
site.site.
(a) Meanwind
speed
in eachwind
direction;
Figure7. Validationof windspeed
at themast(a) Meanwind
speed
in eachwind
direction;(b) Wind speed frequency distribution, (c) Frequency distribution of wind direction.
(b) Wind speed frequency distribution, (c) Frequency distribution of wind direction.

Finally, the proposed microscale flow model overall is validated for all the wind directions. Figure 8 shows comparison of measured and simulated wind speed and winddirection at the met mast location at the hub height. As depicted in the figures and thefitted lines indicate, the microscale flow model results can represent the measurement
well. The mean error of the predicted wind speed at the met mast is −4.7%, which satisfies

Atmosphere 2024, 15, 723

14 of 22

Atmosphere 2024, 15, x FOR PEER REVIEW Finally, the proposed microscale flow model overall is validated for all the wind 15 of 24

directions. Figure 8 shows comparison of measured and simulated wind speed and winddirection at the met mast location at the hub height. As depicted in the figures and thefitted lines indicate, the microscale flow model results can represent the measurement well.
Atmosphere 2024, 15, x FOR PEER REVIEW
15 of 24
The mean error of the predicted wind speed at the met mast is −4.7%, which satisfies theminimum requirement of the mean error of 5%.

(a)

(b)

(a)
(b) and (b) wind direction at the met
Figure 8. Comparison of measured and simulated (a) wind speed
mast hubheight. of measured and simulated (a) wind speed and (b) wind direction at the met
Figure8. Comparison
Figure 8. Comparison of measured and simulated (a) wind speed and (b) wind direction at the met
mast hub height.
mast hub height.

3.3. Modeling of Atmospheric Stability and Its Effect on the Vertical Extrapolation

3.3.
ofofAtmospheric
Stability
ItsItsEffect
Extrapolation
3.3.Modeling
Modeling
Atmospheric
Stabilityand
andcalculated
Effectononthe
theVertical
Vertical
Extrapolation
The pressure
and temperaturebyusing
mesoscalemodel WRF are firstlyThe
pressureandtemperature
calculatedbyusing
mesoscale
model
WRF
The
pressureandtemperature
calculatedbyusing
mesoscalemodel
WRFarearefirstly
firstlyvalidated at Mutsu meteorological station. Figure 9 shows the comparisonof the simuvalidated
at
Mutsu
meteorologicalstation.
Figure9shows
thecomparisonofthesimuvalidated
at
Mutsu
meteorologicalstation.
Figure9shows
thecomparison
ofthesimulated
lated and measured sea level pressure at Mutsu meteorological station. Figure 10a prelated
measuredsea levelpressure
at Mutsu
meteorological station.
Figurepreand andmeasuredsea levelpressure
at Mutsu
meteorological
10a10apresents
sents
the averagetemperatureatthe samehour of the station.
day forFigure
the wholemeasurement
sents
the
averagetemperature
at
the
same
hour
of
the
day
for
the
whole
measurement
the average temperature at the same hour of the day for the whole measurement period.
period.Figure
Figure10b illustrates
the scatter
plot
of1-h
thevalue.
1-h value.
The predictedsurface temperiod.
10b illustrates
the scatterof theThe predicted
temFigure 10b illustrates
the scatterplot ofplot
the 1-hvalue.
The predictedsurfacesurface
temperature
perature
and
sea
level
pressureshow
agreementwiththe measurements.
This result
peraturesea
levelpressureshow
goodgood
agreementthemeasurements.
Thisindicates
resultand sea andlevel
pressureshow
goodagreementwiththewithmeasurements.
This result
indicatesthatthemesoscalesimulationisreliable.
that the thatmesoscalesimulationis reliable.
indicates
the mesoscale
simulationis reliable.

Sea Level Pressure [hPa]

Sea Level Pressure [hPa]

1040
1040
1030
1030
1020

1020

1010

1010

1000

1000

Simulation
Measurement
Simulation

990

990
980
11/23

980
11/23

Measurement
11/28

12/3

12/8

11/28

12/3

Time 12/8(day)

Time (day)

(a)

12/13

12/13

12/18

12/23

12/18

12/23

(b)

(a) 9. Comparison of predicted and measured sea level pressure at the(b) Mutsu meteorologicalFigurestation.
(a)Timeseries of ofhourlydata 23andNovember–23
2018; (b) Scatterof hourly
Figure9. Comparison
predicted
measured December
sea level pressure
at the plot
Mutsu
meteorologicalFigureof predicted
data
from9.20Comparison
November 2018–30
June and2019.measured sea level pressure at the Mutsu meteorologicalstation.
(a)Timeseries
ofhourlydata
23
November–23
December
2018;(b)
Scatterplotof hourly
station. (a) Timeseries of hourly data 23 November–23 December 2018; (b) Scatter plot of hourly
data
data
20 November
2018–30
June 2019.
fromfrom
20 November
2018–30
June 2019.

Atmosphere 2024, 15, x FOR PEER REVIEW

16 of 24

Atmosphere 2024, 15, 723

15 of 22

10
Simulation
Measurement

Temperature [癈]

8

6

4

2

0
0

6

12

18

24

Hour of the day

(a)

(b)

Figure 10. Comparison of predicted and measured temperature at the Mutsu meteorological station.
(a) Average temperature at the same hour of the day for the whole measurement period; (b) Scatterplot of hourly data from 20 November 2018–30 June 2019.

The directionalvariation of atmospheric stability is then investigated using(a)(b) 1/L values extracted from mesoscale simulation. The variation of atmospheric stability with the
10. Comparison
predicted
and temperature
measuredtemperature
the Mutsu
wind
directionis shown
inofFigure
11.ofThis
classificationis conductedaccording
toattheFigure10. Figure
Comparison
predicted
andmeasuredatthe Mutsu
meteorologicalstation.
classes shownin Table 5. Thereis a clear
tendency
towardstable atmosphericconditions,meteorologicalstation.
(a)
Averagetemperature
at
the
same
hour
of
theperiod;day (b) Scatter(a) Average temperature at the same hour of the day for the whole measurement
especially from
S towholeSW wind
where the windis blowing
from plot
landofcompared
to other
for themeasurement
period;(b) Scatter
hourly datafrom 20plot ofhourly datafrom
20 November
2018–30
June
2019.
wind
directions.
This
clear
correlation
betweenwind
directionandatmosphericstability
November 2018–30 June 2019.
implies the justification of the use of directional stability correction factors as defined in
The(37).
directional variation of atmospheric stability is then investigated using 1/L values
Equationextracted from mesoscale simulation. The variation of atmospheric stability with thewind direction is shown in Figure 11. This classification is conducted according to the
Atmospheric
Stability
Monin–Obukhovclasses shownin Table
5. There is a clear tendency
towardLength
stable atmospheric conditions,Very< L < 50 from land compared to otherespecially fromS stableto SW wind where the wind is10blowing
50 < L < 200wind directions.Stable
This clear correlation between winddirection and atmospheric stability
Neutral|L| > 200
implies the justification
of
the
use
of
directional
stability
correction factors as defined in
Unstable−200 < L < −100
Equation (37).
Very unstable−100 < L < −50

Table 5. Classification of atmospheric stability according to Monin–Obukhov length scale.

FigureDirectionalvariabilityof atmosphericstability.
Figure11.11.
Directionalvariabilityof atmosphericstability.

Table 5. Classification of atmospheric stability according to Monin–Obukhov length scale.
Atmospheric Stability

Monin–Obukhov Length

Very stable
StableNeutralUnstable
Very unstable

10 < L < 50
50 < L < 200|L| > 200
−200 < L < −100
−100 < L < −50

The equivalent value of 1/Leq is calculated for each wind direction using
Equation (39) and shown in Figure 12. The low wind speed data does not contributeto the mean wind speed profile and is excluded from this analysis. Three different criteria
are used to filter low wind speed data, i.e., wind speeds lower than 2 m/s [40], 3 m/s [14]1Atmosphere 2024, 15, 723

16 of 22and 4 m/s [39] are filtered. Figure 12 shows the equivalent value of 1/Leq by using thesecriteria. The harmonic average of L proposed by Duran et al. [14] is also plotted in Figure 12.

Figure 12. Directional value of 1/Leq obtained using the proposed method and harmonic average [14].

It is found that the proposed method is not sensitive to the wind speed filter criteria.
The results are an order of magnitude different with the harmonic average of L proposedby Duran et al. [14]. In this study, a 4 m/s filter is used. The vertical profile of mean windspeed in NNW direction at the lidar site is plotted in Figure 13, In this direction, the mean
wind speed is not affected by the turbine wake. The harmonic average of L overestimatesthe effect of stability, while the proposed method shows favorable agreement with themeasurement.

Figure 13. Comparison of the predicted and measured vertical profile of mean wind speed at thelidar site in NNW wind direction [14].

The predicted and measured vertical profiles of mean and standard deviation of thewind speed for sixteen wind directions at the lidar site are shown in Figures 14 and 15.
The hub height as well as the upper and the lower tip of the turbine blade are marked
using dashed lines. It is clear that the predicted vertical profile shows favorable agreementwith the measurement when topography effect, atmospheric stability and turbine wake areconsidered, except for easterly wind where the wind is affected by the wake of the windturbines located on a different ridge, the elevation of which is slightly higher than the lidar
site. Wake effect of such case is difficult to express by the proposed approach. However, thefrequency of occurrence for easterly wind is low and the effect on the overall performanceis minor.

Atmosphere 2024, 15, x FOR PEER REVIEW
Atmosphere 2024, 15, 723

19 of 24
17 of 222.31 %

0.57 %

0.44 %

0.26 %

0.57 %

1.43 %

5.22 %

3.14 %

1.65 %

2.42 %

5.23 %

9.56 %

21.4 %

29.47 %

13.26 %

3.06 %

Comparison ofof predictedpredicted andand measured
measured vertical profiles of mean wind speed at the lidar
Figure 14. Comparisonsite for 16 wind directions.

When the effect of atmospheric stability is considered (blue line), the overestimation ispartly improved (e.g., SSE to WSW). Consideration of the turbine wake further improves
the predicted mean wind speeds and turbulence intensities (green line), especially in thewind directions when the lidar is in some of the turbines’ wakes (e.g., NNE, NE, E, ESE,
SE, SSE, S, SSW and SW). The overestimation of the wind speed and underestimation of
turbulence intensity are significantly improved by considering the effect of turbine wake.

Atmosphere 2024, 15, 723

18 of 22

Figure 15. Comparison of predicted and measured vertical profile of turbulence intensity at the lidarsite for 16 wind directions.

4. Conclusions
In this study, a microscale flow model is proposed considering the effect of topography,
atmospheric stability, and wind turbine wake, and validated by using measurement data at
the Iwaya Wind Farm, north of Japan. The following conclusions are obtained.
1.

The predicted wind speed and direction at the met mast by using lidar-measured
wind velocity and the proposed microscale flow model shows good agreement withthe cup anemometer measurement at the met mast when the effect of topography,
atmospheric stability and wind turbine wake are considered.

Atmosphere 2024, 15, x FOR PEER REVIEW

Atmosphere 2024, 15, 723

21 of 24

19 of 22the cup anemometer measurement at the met mast when the effect of topography,
atmospheric stability and wind turbine wake are considered.
2. A stability correction factor is proposed using the Monin–Obukhov length scale ob2.
A stability correction factor is proposed using the Monin–Obukhov length scale
tained from the mesoscale model to evaluate the effect of atmospheric stability on theobtained from the mesoscale model to evaluate the effect of atmospheric stability onmean wind speed. The proposed equivalent Monin–Obukhov length scale gives betthe mean wind speed. The proposed equivalent Monin–Obukhov length scale gives
ter agreement with the measurements than the conventional harmonic average ofbetter agreement with the measurements than the conventional harmonic average ofMonin–Obukhov length scale.
Monin–Obukhov length scale.
3. The effects of topography, atmospheric stability, and turbine wake on the vertical proThe effects of topography, atmospheric stability, and turbine wake on the vertical3.
file of wind speed and turbulence intensity are also investigated. The consideration of
profile of wind speed and turbulence intensity are also investigated. The consideration
all threegives betterwith themeasurements
by theby
verticallidar. lidar.
of alleffects
three effects
gives agreement
better agreementwiththe measurements
the vertical
Author Contributions: Conceptualization, A.Y., A.T. and T.I.; methodology, A.Y., A.T. and T.I.;
AuthorContributions:Conceptualization,
and T.I.; methodology,
A.T. and T.I.;
software,A.T.;
validation, A.T.;
formal analysis, A.Y.,
A.T.; A.T.
investigation,
A.Y., A.T. and A.Y.,
T.I.; resources,
validation,
A.T.;
formalanalysis,
A.T.; investigation,
A.Y., A.T. andA.Y.software,
and T.I.; A.T.;
data curation,
A.Y.,
A.T.
and T.I.;
writing—originaldraft preparation,
A.Y.T.I.;
andresources,
A.T.;
A.Y. and T.I.; anddataediting,curation,
A.Y.,
A.T.
and
writing—originaldraft preparation,
writing—review
A.Y.,
A.T.
and
T.I.;T.I.;
visualization,
A.T.; supervision,
A.Y. andA.Y.
T.I.;andpro-A.T.;
andandediting,
A.Y., A.T.
and T.I.; visualization,
supervision,
A.Y. andproject
ject writing—review
administration A.Y.
T.I.; funding
acquisition,
T.I. and A.Y.A.T.;
All authors
have readand T.I.;
agreedadministration
A.Y. andT.I.;manuscript.
funding acquisition, T.I. and A.Y. All authors have read and agreed toto thepublished versionof thethe published version of the manuscript.
Funding: This research received no external funding.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: The data presented in this study are available on request from theData Availability
presented
in thisstudy
areavailable
on requestfrom thecorresponding
author. Statement:
The data areThenotdatapublicly
available
becausetheauthors
do nothave permissioncorresponding
to share data. author. The data are not publicly available because the authors do not have permission
to share data.
Acknowledgments: The authors would like to acknowledge Eurus Energy Holdings Corporation
for themet mast and lidar
at thewould
Iwayalike
Wind
Farm.
Acknowledgments:
Thedataauthors
to acknowledge
Eurus Energy Holdings Corporation forthe met mast and lidar data at the Iwaya Wind Farm.
Conflicts of Interest: The authors declare no conflicts of interest.
Conflicts of Interest: The authors declare no conflicts of interest.

Appendix A. Correction of Lidar-Measured Wind Speed over Complex Terrain
Appendix
A. Correctionof Lidar-Measured
Wind 𝑉SpeedoverComplexIn the pulsed
lidar, the line-of-sightwind speed, at anyheighthas aTerrain
geometrical𝑖

In thewith
pulsedline-of-sight
any height
has a geometrical
relationshipthe lidar,
windthespeed
vector 𝐮wind) i ,asat shownin FigureA1 and𝑖 = (𝑢speed
𝑖 , 𝑣𝑖 , 𝑤𝑖 V
relationship
Equation(A1). with the wind speed vector ui = (ui , vi , wi ) as shown in Figure A1 andEquation (A1).
(A1) (A1)
𝑉𝑖 = 𝐮V𝑖 ∙=𝐧u= ·𝑢n𝑖 cos𝜑sin𝜃
𝑣𝑖 sin𝜑sin𝜃
𝑖 cos𝜃
= ui cos φ+sin
θ + vi sin+φ 𝑤
sin
θ + wi cos θiiFigureA1. A1.
Schematicof themeasurement.
FigureSchematicof lidarthe lidar
measurement.

Assuming
thatthat
the thewindfieldfield
is homogenous
on theplaneat anyheight,the thewind
Assumingwind
is homogenous
on theplaneat anyheight,wind
speed
accordingto theDoppler
Beam
Swing
(DBS)
methodusing
the theline-of-sight
wind
speed
accordingto theDoppler
Beam
Swing
(DBS)
methodusing
line-of-sightwind
speedscanbecalculatedbyEquations
(A2)
and
(A3).
speeds can be calculated by Equations (A2) and (A3).
𝑢dbs =u

𝑉2 − 𝑉4V − V
4
=𝜃 2
dbs
2 sin
2 sin θ(A2) (A2)V1 − V3
2 sin θ(A3)

vdbs =

Atmosphere 2024, 15, x FOR PEER REVIEW

22 of 24

Atmosphere 2024, 15, 723

𝑣dbs =

𝑉1 − 𝑉3
2 sin 𝜃20 of 22

(A3)

However, this assumption on complex terrain is not correct, and the measured wind
However, this assumption on complex terrain is not correct, and the measured windspeed is erroneous. CFD were used to calculate the wind speed according to the DBS
speed is erroneous. CFD were used to calculate the wind speed according to the DBS
method as shown in [27].
method as shown in [27].
The horizontal wind speed measured by the lidar is calculated as follows:
The horizontal wind speed measured by the lidar is calculated as follows:q
𝑉dbs = √𝑢2dbs + 𝑣2dbs
Vdbs = u2dbs + v2dbs

(A4)

(A4)

Then, the correction factor 𝛾 can be calculated as follows:
Then, the correction factor γ can be calculated as follows:(A5)
𝛾 = 𝑉𝑐 ⁄𝑉dbs
γ = Vc /Vdbs(A5)
where 𝑉𝑐 is the value of the wind speed in the middle of the circle in the CFD simulations.
FigureA2 shows
thevaluecorrection
threeheights
as acircle
functionof windwhere
Vc is theof thefactors
wind forspeed
in typicalthe middleof thein the
CFD disimulations. Figure A2 shows the correction factors for three typical heights as a function ofrection.
wind
direction.
Finally,
the predicted and measured wind speed at the met mast hub height using
Finally,
predicted
andaftermeasuredwind
at the metfactor
mast hubheightinusing
lidar 160-m
wind thespeed
before andapplying
thespeedlidar correction
are shown
lidar
speedthebeforeapplying the lidar
correction
factorarewindshown inFigureA3.160-m
In thiswind
example,slopeandandafterthe determination
coefficientfor themean
FigureA3. In thisand thedetermination coefficient for the mean wind
speed
are improved
byexample,
applyingthe
theslope
correction
factor.
speed are improved by applying the correction factor.
1.1

Correction factor

z=40.3 m
z=100 mz=159.8 m
1.05

1

0.95
0

60

120

180

240

300

360

Wind Direction [degree]
Atmosphere 2024, 15, x FOR PEER REVIEW

23 of 24

Figure A2. Correction factor for the lidar-measured wind speed at three typical heights at the
Figure A2. Correction factor for the lidar-measured wind speed at three typical heights at the lidar
lidar location.
location.

(a)

(b)

Figureof
A3.
Comparison
simulated andmeasuredspeed
at the(a)metmast. applying(a) Before applying
Figure A3. Comparison
simulated
andofmeasured
wind
speed atwind
the metmast.
Beforethe lidar correction factor; (b) After applying the lidar correction factor.
the lidar correction factor; (b) After applying the lidar correction factor.

References
1.
2.
3.

Li, S.; Sun, X.; Zhang, S.; Zhao, S.; Zhang, R. A study on microscale wind simulations with a coupled WRF–CFD model in the
Chongli mountain region of Hebei Province, China. Atmosphere 2019, 10, 731. https://doi.org/10.3390/atmos10120731.
Tomaszewski1, J.M.; Lundquist, J.K. Simulated wind farm wake sensitivity to configuration choices in the Weather Researchand Forecasting model version 3.8.1. Geosci. Model Dev. 2020, 13, 2645–2662. https://doi.org/10.5194/gmd-13-2645-2020.
WindFarmer, G.H. Wind Farm Design Software 4.0: Theory Manual; Garrad Hassan and Partners Ltd.: Bristol, UK, 2010.

Atmosphere 2024, 15, 723

21 of 22

References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.

21.
22.
23.
24.
25.
26.
27.
28.
29.
30.

Li, S.; Sun, X.; Zhang, S.; Zhao, S.; Zhang, R. A study on microscale wind simulations with a coupled WRF–CFD model in the
Chongli mountain region of Hebei Province, China. Atmosphere 2019, 10, 731. [CrossRef]
Tomaszewski1, J.M.; Lundquist, J.K. Simulated wind farm wake sensitivity to configuration choices in the Weather Research andForecasting model version 3.8.1. Geosci. Model Dev. 2020, 13, 2645–2662. [CrossRef]
GH WindFarmer. Wind Farm Design Software 4.0: Theory Manual; Garrad Hassan and Partners Ltd.: Bristol, UK, 2010.
IEC 61400-12-1; Wind Turbines—Part 12-1: Power Performance Measurements of Electricity Producing Wind Turbines, Edition
1.0. International Electromechanical Commission: Geneva, Switzerland, 16 December 2005.
Murthy, K.S.R.; Rahi, O.P. A Comprehensive Review of Wind Resource Assessment. Renew. Sustain. Energy Rev. 2017, 72,
1320–1342. [CrossRef]
Ross, A.N.; Arnold, S.; Vosper, S.B.; Mobbs, S.D.; Dixon, S.D.; Robins, A.G. A Comparison of Wind-Tunnel Experiments andNumerical Simulations of Neutral and Stratified Flow Over a Hill. Bound.-Layer Meteorol. 2004, 11, 427–459. [CrossRef]
Bechmann, A.; Sørensen, N.N. Hybrid RANS/LES applied to complex terrain. Wind Energy 2011, 14, 225–237. [CrossRef]
Bleeg, J.; Digraskar, D.; Woodcock, J.; Corbett, J.F. Modeling stable thermal stratification and its impact on wind flow overtopography. Wind Energy 2015, 18, 369–383. [CrossRef]
Meissner, C.; Gravdahl, A.R.; Steensen, B. Including Thermal Effects in CFD Wind Flow Simulations. J. Environ. Sci. Int. 2009, 18,
833–839. [CrossRef]
Liu, L.; Stevens, R.J. Effects of atmospheric stability on the performance of a wind turbine located behind a three-dimensional hill.
Renew. Energy 2021, 175, 926–935. [CrossRef]
Qian, G.W.; Ishihara, T. A Novel Probabilistic Power Curve Model to Predict the Power Production and Its Uncertainty for aWind Farm over Complex Terrain. Energy 2022, 261, 125171. [CrossRef]
Corbett, J.-F.; Poenario, A.; Horn, U.; Leask, P. An Extensive Validation of CFD Flow Modelling. In Proceedings of the DEWEK,
Bremen, Germany, 18–19 May 2015.
Uchida, T.; Takakuwa, S. Numerical investigation of stable stratification effects on wind resource assessment in complex terrain.
Energies 2020, 13, 6638. [CrossRef]
Durán, P.; Meiβner, C.; Casso, P. A New Meso-Microscale Coupled Modelling Framework for Wind Resource Assessment: A
Validation Study. Renew Energy 2020, 160, 538–554. [CrossRef]
Stiperski, I.; Calaf, M. Generalizing Monin-Obukhov Similarity Theory (1954) for Complex Atmospheric Turbulence. Phys. Rev.
Lett. 2023, 130, 124001. [CrossRef] [PubMed]
Qian, G.-W.; Ishihara, T. Wind Farm Power Maximization through Wake Steering with a New Multiple Wake Model for Predictionof Turbulence Intensity. Energy 2021, 220, 119680. [CrossRef]
Ishihara, T.; Qian, G.W.; Qi, Y.H. Numerical Study of Turbulent Flow Fields in Urban Areas by Using Modified k-ε Model andLarge Eddy Simulation. J. Wind. Eng. Ind. Aerodyn. 2020, 206, 1–20. [CrossRef]
Oke, T.R. Boundary Layer Climates; Routledge: London, UK, 2002; ISBN 1134951337.
Ferziger, J.; Peric, M. Computational Method for Fluid Dynamics, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2002.
Skamarock, W.C.; Klemp, J.B.; Dudhia, J.B.; Gill, D.O.; Barker, D.M.; Duda, M.G.; Huang, X.-Y.; Wang, W.; Powers, J.G. A
Description of the Advanced Research WRF Version 3, NCAR Technical Note TN-475+STR. 2008. Available online: https://opensky.ucar.edu/islandora/object/technotes:500 (accessed on 10 November 2023).
Stauffer, D.; Seaman, N.L. Use of Four-Dimensional Data Assimilation in a Limited-Area Mesoscale Model. Part I: Experimentswith Synoptic-Scale Data. Mon. Weather Rev. 1990, 118, 1250–1277. [CrossRef]
Webb, E.K. Profile Relationships: The Log-Linear Range, and Extension to Strong Stability. Q. J. R. Meteorol. Soc. 1970, 96, 67–90.
[CrossRef]
Dyer, A.J. A Review of Flux-Profile Relationships. Bound. Layer Meteorol. 1974, 7, 363–372. [CrossRef]
Holtslag, A.A.M.; De Bruin, H.A.R. Applied Modeling of the Nighttime Surface Energy Balance over Land. J. Appl. Meteorol. 1988,
27, 689–704. [CrossRef]
Paulson, C.A. The Mathematical Representation of Wind Speed and Temperature Profiles in the Unstable Atmospheric Surface
Layer. J. Appl. Meteorol. 1970, 9, 857–861. [CrossRef]
Klaas-Witt, T.; Emeis, S. The Five Main Influencing Factors for Lidar Errors in Complex Terrain. Wind. Energy Sci. 2022, 7, 413–431.
[CrossRef]
Vogstad, K.; Simonsen, A.H.; Brennan, K.J.; Lund, J.A. Uncertainty of Lidars in Complex Terrain. In Proceedings of the European
Wind Energy Conference and Exhibition, Vienna, Austria, 4–7 February 2013.
Ishihara, T.; Yamaguchi, A.; Fujino, Y. A nonlinear model for predictions of turbulent flow over steep terrain. In Proceedings ofthe First World Wind Energy Conference and Exhibition, VB3.4, Berlin, Germany, 4–8 July 2002; pp. 1–4.
Qian, G.W.; Ishihara, T. Numerical study of wind turbine wakes over escarpments by a modified delayed detached eddy
simulation. J. Wind. Eng. Ind. Aerodyn. 2019, 191, 41–53. [CrossRef]
Stark, J.D.; Donlon, C.J.; Martin, M.J.; McCulloch, M.E. OSTIA: An Operational, High Resolution, Real Time, Global Sea Surface
Temperature Analysis System. In Proceedings of the OCEANS 2007—Europe, Aberdeen, UK, 18–21 June 2007; pp. 1–4. [CrossRef]

Atmosphere 2024, 15, 723

31.

32.

33.
34.
35.

36.
37.
38.
39.

40.

22 of 22

Hersbach, H.; Bell, B.; Berrisford, P.; Biavati, G.; Horányi, A.; Muñoz Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Rozum, I.; et al.
ERA5 Hourly Data on Single Levels from 1940 to Present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS).
2023. Available online: https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24381/cds.adbb2d47?tab=overview (accessed on
12 June 2024).
Ferrier, B.S.; Jin, Y.; Lin, Y.; Black, T.; Rogers, E.; Dimego, G. Implementation of a New Grid-Scale Cloud and Precipitation
Scheme in the NCEP Eta Model. In Proceedings of the 19th Conference on weather Analysis and Forecasting/15th Conference on
Numerical Weather Prediction, San Antonio, TX, USA, 12–16 August 2002; pp. 280–283.
Iacono, M.J.; Delamere, J.S.; Mlawer, E.J.; Shephard, M.W.; Clough, S.A.; Collins, W.D. Radiative Forcing by Long-Lived
Greenhouse Gases: Calculations with the AER Radiative Transfer Models. J. Geophys. Res. 2008, 113, D13103. [CrossRef]
Dudhia, J. Numerical Study of Convection Observed during the Winter Monsoon Experiment Using a Mesoscale Two-Dimensional
Model. J. Atmos. Sci. 1989, 46, 3077–3107. [CrossRef]
Tewari, M.; Chen, F.; Wang, W.; Dudhia, J.; Lemone, M.A.; Mitchell, K.; Ek, M.; Gayno, G.; Wegiel, J.; Cuenca, R.H. Implementation
and Verification of the United NOAH Land Surface Model in the WRF Model. In Proceedings of the 20th Conference on Weather
Analysis and Forecasting/16th Conference on Numerical Weather Prediction, Seattle, WA, USA, 12–15 January 2004; pp. 11–15.
Mellor, G.L.; Yamada, T. Development of a Turbulence Closure Model for Geophysical Fluid Problems. Rev. Geophys. 1982,
20, 851. [CrossRef]
Betts, A.K.; Miller, M.J. A New Convective Adjustment Scheme. Part II: Single Column Tests Using GATE Wave, BOMEX, ATEX
and Arctic Air-Mass Data Sets. Quart. J. Roy. Meteor. Soc. 1986, 112, 693–709. [CrossRef]
Kikuchi, Y.; Fukushima, M.; Ishihara, T. Assessment of a Coastal Offshore Wind Climate by Means of Mesoscale Model Simulations
Considering High-Resolution Land Use and Sea Surface Temperature Data Sets. Atmosphere 2020, 11, 379. [CrossRef]
The Carbon Trust. OWA Roadmap for the Commercial Acceptance of Floating LiDAR Technology. 2018. Available online: https://www.carbontrust.com/our-work-and-impact/guides-reports-and-tools/roadmap-for-commercial-acceptanceof-floating-lidar (accessed on 1 September 2023).
IEA Wind. Expert Group Study on Recommended Practices, 15. Ground-Based Vertically-Profiling Remote Sensing for Wind
Resource Assessment, First Edition. 2013. Available online: https://www.fulcrum3d.com/wp-content/uploads/2018/01/
Ground-Based-Vertically-Profiling-Remote-Sensing-for-Wind-Assessment.pdf (accessed on 1 September 2023).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.