Assessment of Wind over Complex Terrain Considering the Effects of Topography, Atmospheric Stability and Turbine Wakes

: This study proposes a microscale flow model to estimate mean wind speed, fluctuating wind speed and wind direction over complex terrain considering the effects of topography, atmospheric stability, 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 atmospheric stability. 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 to account for the effect of atmospheric stability. Finally, the measured wind data at the reference height are vertically extrapolated and validated at the lidar location. The predicted mean and fluctuating wind speeds show good agreement with the measurements.


Introduction
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 wind prediction 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 wind turbine positions.According to the WindFarmer theory manual [3], three main sources of 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 wind condition 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 of wind 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][7][8][9][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 the temperature field.However, there are still a few problems to be solved.
The first problem is the cross validation of the flow model for the repowering of the existing wind farms and development of new wind farms near the existing ones.In this situation, almost all wind direction is somehow affected by the wake of the existing wind turbines and this wake effect must be included in the flow model.In addition, in the site where the atmospheric stability is important, this effect also must be included.When the surface 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 cross validation 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 computational cost, 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 of Monin-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 wake model on the vertical profile of mean wind speed and turbulence intensity need to be considered 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 wake model is introduced to account for the effects of turbine wakes.The predicted mean and fluctuating 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.

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 by considering 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 caused by the topography and surface roughness.The details of the CFD model are described in Section 2.1.If measured mean wind speed, U obs ref , fluctuating wind speed, σ obs ref , and wind direction, θ obs ref are available at one point, by using terrain correction factors, these variables at any point (x, y, z) can be calculated.
where C terrain U , C terrain σ are the terrain correction factors for the mean and fluctuating wind speeds, respectively, and ∆θ terrain is the correction term for the wind direction and x = (x, y, z) is the position of the site where the wind is calculated.
where   t ,   t are the terrain correction factors for the mean and fluctuating wind speeds, respectively, and ∆ t is the correction term for the wind direction and  = (, , ) is the position of the site where the wind is calculated.Then, the Monin-Obukhov length scale  is extracted from the mesoscale simulation and used to calculate the stability correction factor   stability (,  ref obs ) which is also a function of the point () and wind direction  ref obs .The detail of this correction factor is described in Section 2.2.According to Stiperski and Calaf [15], the effect of stability on fluctuating wind speed is small in neutral, stable, and moderately unstable cases.So, in this study, the effect of stability on the fluctuating wind speed is neglected.

Onsite measurement
Terrain correction factor for mean wind speed  (, ) Terrain correction factor for fluctuating wind speed  (, ) Terrain correction term for wind direction  (, )

CFD simulation
Wake free mean wind speed  target initial (, ) Wake free fluctuating wind speed  target initial (, ) Wake free wind direction  target initial (, )

Mesoscale model
Terrain corrected mean wind speed  target terrain (, ) Then, the Monin-Obukhov length scale L is extracted from the mesoscale simulation and used to calculate the stability correction factor C stability U x, θ obs ref which is also a function of the point (x) and wind direction θ obs ref .
The detail of this correction factor is described in Section 2.2.According to Stiperski and Calaf [15], the effect of stability on fluctuating wind speed is small in neutral, stable, and moderately unstable cases.So, in this study, the effect of stability on the fluctuating wind speed is neglected.
Finally, an advanced wake model is applied to calculate the velocity deficit ∆U(x, y, z, t), and added turbulence ∆σ(x, y, z, t) from wind turbines.The detail of the wake model is described in Section 2.3.
where U pred target (x, t), I pred target (x, t) and θ pred target (x, t) are the predicted mean wind speed, tur- bulence intensity and wind direction, respectively.The effect of multiple turbines is superposed 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 the mean wind speed can be obtained.
This equation is based on dimensional analysis and the effects of terrain and atmospheric stability are accounted for by multiplying the non-dimensional coefficients C terrain .This paper aims to propose a method to identify these coefficients.

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.
where u i is the filtered wind speed in the x i direction.ρ is the air density, p is the filtered pressure and µ is the molecular viscosity of the air.τ ij is the SGS stresses which is modeled as follows: where µ t denotes the SGS turbulent viscosity, S ij 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: where L s 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 C s , 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).
where c f is a kind of drag coefficient representative, l 0 is the representative length scale of the obstacle and γ 0 is the packing density.The values of c f , l 0 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 is shifted 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 the second 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 wind speed, and wind direction at any point inside the domain using reference point data by correcting the effect of terrain and roughness in each direction sector.The correction factor for each of these variables is defined in Equations ( 18)- (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 direction sector.The correction factors for any direction, θ obs ref , are obtained using a linear interpolation.
First, the direction sector of θ obs ref is obtained as shown in Equation (21).
The linear interpolation is then performed as follows.
where the proportional coefficient a is defined in Equation (24).

Mesoscale Model and Stability Correction
As a mesoscale model, the Weather Research and Forecasting (WRF, 3.4.1)modeling system, 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 hundreds of 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 the equations are listed in Table 1.

Types Parameter Definition
Prognostic variables Moist potential temperature q v , q c , q r , q i mixing ratios for water vapor, cloud, rain, ice Diagnostic variables Hydrostatic component of pressure for dry air P ht P h along the top boundary for dry air The ARW equations are formulated using a terrain-following hydrostatic pressure vertical coordinate denoted by η.
where µ d is the mass of dry air in the column and P dh and P dht are hydrostatic pressure of the dry 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: The governing thermodynamics equation, conservation of mass, geopotential time derivative and the equation of state are, respectively, shown in Equations ( 30)- (33).
where α is the inverse density taking the whole parcel into account 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 as shown in Equation ( 35) [21].
where θ is the prognostic variable, F(θ) represents the normal tendency terms due to physics, the nudging parameter, G θ is the timescale controlling the nudging strength and θ 0 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 atmospheric stability 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][23][24][25][26].
where u * is the friction velocity, κ is the von Karman constant, z 0 is the surface roughness and L is the Monin-Obukhov length scale, which is extracted from the mesoscale simulation results; the function ψ z L is defined in Equation (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: The idea for Equation (37) is that the effect of topography and the effect of atmospheric stability can be separated, and the stability correction factor C stability U_t (x, t) is defined assuming the terrain is flat.The value C stability U_t (x, t) becomes unity for the same height (z = z ref ) 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 an equivalent 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 is proposed to calculate the directional stability correction factor.The proposed factor is calculated using Equation (39).(39) where ψ eq (x, n θ ) can be calculated as shown in Equation (39).
U obs ref (t) is the observed wind speed in that direction sector at time t.Once ψ eq is calculated for each direction sector, the equivalent value of 1/L eq for this wind direction can be calculated inversely which can explain the effect of stability for this wind direction and expressed as: To obtain C stability U x, θ obs ref at any direction θ obs ref , the C stability U (x, n θ ) is linearly interpo- lated in the similar way as explained in Section 2.1.

Turbine Wake Model
In 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 that has 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.

Wake Model Formula Parameters
Wake width where D is the rotor diameter, C T is the thrust coefficient, U h is the rotor onset wind speed and I a is the rotor onset turbulence intensity.
As the turbines experience a non-uniform flow caused by the complex terrain or the wake of wind turbines further upstream, it is first required to calculate the rotor onset wind speed U h and turbulence I a .The rotor onset wind speed U h is calculated by directly performing a geometric averaging of wind speed U over the rotor, while rotor onset I a is calculated by the root mean of squares of streamwise turbulence standard deviation σ u over the rotor divided by the rotor onset wind speed U h [16].
Using these values, the effect of wake on the flow field, ∆U and ∆σ, can be calculated for each wind turbine.

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.

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 At the Met mast and lidar locations, several wind direction sectors are affected by turbine wakes.Figure 3 shows the alignment of turbines relative to the lidar and met mast.The wind rose and frequency distribution of the wind speed at lidar 160 m are shown in Figure 4.
The computation domain of the CFD simulation is shown in Figure 5a and the mesh is 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 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 explained in 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 size of 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  mesh size of 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 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 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 defined so 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 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 ratio of 1.1 and the top boundary is located 10,000 m above sea level.The computation domain of the CFD simulation is shown in Figure 5a and the mesh is 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 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 explained in 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 size of 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 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 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 defined so 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 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 ratio of 1.1 and the top boundary is located 10,000 m above sea level.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 and buffer zones for the boundary conditions are added making the total computational domain 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 ratio of 1.1 and the top boundary is located 10,000 m above sea level.Table 3. CFD simulation domain and mesh setting.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.

Domain and Mesh Setting
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.

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) and wind direction (70 m above ground level) are used as reference data and wind speed and direction at the met mast are predicted by using the proposed method.Very low wind speeds cause large uncertainties and are not important in power generation.Therefore, wind speeds lower than 4 m/s [39] are filtered out.
The predicted mean wind speed at the met mast at the hub height is plotted as a function of wind direction and compared with the measurement data in Figure 7a.The red line 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 (blue line), 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 mesoscale model 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 is clearly 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 data compared to WRF with the microscale model as shown in Figure 7c.It should be noted that wake model and stability model does not change the wind direction as shown in Equations ( 1)-( 9), and predicted frequency distribution of wind directions are identical for the three models.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 wind direction at the met mast location at the hub height.As depicted in the figures and the fitted 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 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 wind direction at the met mast location at the hub height.As depicted in the figures and the fitted 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 the minimum requirement of the mean error of 5%.

Modeling of Atmospheric Stability and Its Effect on the Vertical Extrapolation
The pressure and temperature calculated by using mesoscale model WRF are firstly validated at Mutsu meteorological station.Figure 9 shows the comparison of the simulated and measured sea level pressure at Mutsu meteorological station.Figure 10a presents the average temperature at the same hour of the day for the whole measurement period.Figure 10b illustrates the scatter plot of the 1-h value.The predicted surface temperature and sea level pressure show good agreement with the measurements.This result indicates that the mesoscale simulation is reliable.

Modeling of Atmospheric Stability and Its Effect on the Vertical Extrapolation
The pressure and temperature calculated by using mesoscale model WRF are firstly validated at Mutsu meteorological station.Figure 9 shows the comparison of the simulated and measured sea level pressure at Mutsu meteorological station.Figure 10a presents the average temperature at the same hour of the day for the whole measurement period.Figure 10b illustrates the scatter plot of the 1-h value.The predicted surface temperature and sea level pressure show good agreement with the measurements.This result indicates that the mesoscale simulation is reliable.

Modeling of Atmospheric Stability and Its Effect on the Vertical Extrapolation
The pressure and temperature calculated by using mesoscale model WRF are firstly validated at Mutsu meteorological station.Figure 9 shows the comparison of the simulated and measured sea level pressure at Mutsu meteorological station.Figure 10a presents the average temperature at the same hour of the day for the whole measurement period.Figure 10b illustrates the scatter plot of the 1-h value.The predicted surface temperature and sea level pressure show good agreement with the measurements.This result indicates that the mesoscale simulation is reliable.The directional variation of atmospheric stability is then investigated using 1/L values extracted from mesoscale simulation.The variation of atmospheric stability with the wind direction is shown in Figure 11.This classification is conducted according to the classes shown in Table 5.There is a clear tendency toward stable atmospheric conditions, especially from S to SW wind where the wind is blowing from land compared to other wind directions.This clear correlation between wind direction and atmospheric stability implies the justification of the use of directional stability correction factors as defined in Equation (37).The directional variation of atmospheric stability is then investigated using 1/L values extracted from mesoscale simulation.The variation of atmospheric stability with the wind direction is shown in Figure 11.This classification is conducted according to the classes shown in Table 5.There is a clear tendency toward stable atmospheric conditions, especially from S to SW wind where the wind is blowing from land compared to other wind directions.This clear correlation between wind direction and atmospheric stability implies the justification of the use of directional stability correction factors as defined in Equation (37).−100 < L < −50 The equivalent value of 1/L eq is calculated for each wind direction using Equation (39) and shown in Figure 12.The low wind speed data does not contribute to 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] Figure 12.Directional value of 1/L eq 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 proposed by Duran et al. [14].In this study, a 4 m/s filter is used.The vertical profile of mean wind speed 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 overestimates the effect of stability, while the proposed method shows favorable agreement with the measurement.The predicted and measured vertical profiles of mean and standard deviation of the wind 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 agreement with the measurement when topography effect, atmospheric stability and turbine wake are considered, except for easterly wind where the wind is affected by the wake of the wind turbines 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, the frequency of occurrence for easterly wind is low and the effect on the overall performance is minor.When the effect of atmospheric stability is considered (blue line), the overestimation is partly 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 the wind 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.

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 with the cup anemometer measurement at the met mast when the effect of topography, atmospheric stability and wind turbine wake are considered.
However, this assumption on complex terrain is not correct, and the measured wind speed is erroneous.CFD were used to calculate the wind speed according to the DBS method as shown in [27].
The horizontal wind speed measured by the lidar is calculated as follows: Then, the correction factor γ can be calculated as follows: where V c is the value of the wind speed in the middle of the circle in the CFD simulations.Figure A2 shows the correction factors for three typical heights as a function of wind direction.Finally, the predicted and measured wind speed at the met mast hub height using lidar 160-m wind speed before and after applying the lidar correction factor are shown in Figure A3.In this example, the slope and the determination coefficient for the mean wind speed are improved by applying the correction factor.
Atmosphere 2024, 15, x FOR PEER REVIEW 22 of 24 dbs =  1 −  3 2 sin  (A3) However, this assumption on complex terrain is not correct, and the wind speed is erroneous.CFD were used to calculate the wind speed according to the DBS method as shown in [27].
The horizontal wind speed measured by the lidar is calculated as follows: dbs = √ dbs 2 +  dbs 2 (A4) Then, the correction factor  can be calculated as follows: where   is the value of the wind speed in the middle of the circle in the CFD simulations.Figure A2 shows the correction factors for three typical heights as a function of wind direction.
Finally, the predicted and measured wind speed at the met mast hub height using lidar 160-m wind speed before and after applying the lidar correction factor are shown in Figure A3.In this example, the slope and the determination coefficient for the mean wind speed are improved by applying the correction factor.

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

Figure 1 .
Figure 1.A framework of the microscale flow model.
These wake-affected directions in which wind turbines exist within 10D, where D = 62 m is the rotor diameter in this wind farm, are shaded in the figure.

Figure 3 .
Figure 3.The placement of the wind turbines relative to (a) Met mast and (b) Lidar.Wake-affected directions are shaded.The wind rose and frequency distribution of the wind speed at lidar 160 m are shown in Figure 4.

Figure 4 .
Figure 4. Wind at Iwaya wind farm lidar.(a) Wind rose; (b) Frequency distribution of the wind speed.

Figure 4 .
Figure 4. Wind at Iwaya wind farm lidar.(a) Wind rose; (b) Frequency distribution of the wind speed.

Figure 5 .
Figure 5. CFD simulation.(a) Domain and boundary conditions; (b) Horizontal mesh on the bottom surface; (c) Typical mesh and the location of the lidar and the met mast for westerly wind direction.

Table 4 .
Description of mesoscale model and computational conditions.

Table 4 .
Description of mesoscale model and computational conditions.

Figure 7 .
Figure 7. Validation of wind speed at the mat mast site.(a) Mean wind speed in each wind direction; (b) Wind speed frequency distribution, (c) Frequency distribution of wind direction.

Figure 7 .
Figure 7. Validation of wind speed at the mat mast site.(a) Mean wind speed in each wind direction; (b) Wind speed frequency distribution, (c) Frequency distribution of wind direction.

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

Figure 9 .
Figure 9.Comparison of predicted and measured sea level pressure at the Mutsu meteorological station.(a) Timeseries of hourly data 23 November-23 December 2018; (b) Scatter plot of hourly data from 20 November 2018-30 June 2019.

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

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

Figure 9 .
Figure 9.Comparison of predicted and measured sea level pressure at the Mutsu meteorological station.(a) Timeseries of hourly data 23 November-23 December 2018; (b) Scatter plot of hourly data from 20 November 2018-30 June 2019.

Figure 9 .Figure 10 .
Figure 9.Comparison of predicted and measured sea level pressure at the Mutsu meteorological station.(a) Timeseries of hourly data 23 November-23 December 2018; (b) Scatter plot of hourly data from 20 November 2018-30 June 2019.

Figure 10 .
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) Scatter plot of hourly data from 20 November 2018-30 June 2019.

Figure 10 .
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) Scatter plot of hourly data from 20 November 2018-30 June 2019.

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

Figure 14 .
Figure 14.Comparison of predicted and measured vertical profiles of mean wind speed at the lidar site for 16 wind directions.

Figure 14 .
Figure 14.Comparison of predicted and measured vertical profiles of mean wind speed at the lidar site for 16 wind directions.

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

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

Figure A2 .
Figure A2.Correction factor for the lidar-measured wind speed at three typical heights at the lidar location.Atmosphere 2024, 15, x FOR PEER REVIEW 23 of 24

Figure A3 .
Figure A3.Comparison of simulated and measured wind speed at the met mast.(a) Before applying the lidar correction factor; (b) After applying the lidar correction factor.

Table 2 .
Summary of the Ishihara-Qian wake model.

Table 3 .
CFD simulation domain and mesh setting.

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

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