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

Abstract

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.

1. 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.

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 by considering the effects of topography, atmospheric stability, and wind turbine wakes as described in Section 2.1, Section 2.2 and Section 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_{\mathrm{ref}}^{\mathrm{obs}}$, fluctuating wind speed, ${\sigma }_{\mathrm{ref}}^{\mathrm{obs}}$, and wind direction, ${\theta }_{\mathrm{ref}}^{\mathrm{obs}}$ are available at one point, by using terrain correction factors, these variables at any point $\left(x,y,z\right)$ can be calculated.

$$U_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)=U_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)\times C_U^{\mathrm{terrain}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$$

(1)

$${\sigma }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)={\sigma }_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)\times C_{\sigma }^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$$

(2)

$${\theta }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)={\theta }_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)+\Delta {\theta }^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$$

(3)

where $C_U^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, $C_{\sigma }^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ are the terrain correction factors for the mean and fluctuating wind speeds, respectively, and $\Delta {\theta }^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ is the correction term for the wind direction and $x=\left(x,y,z\right)$ is the position of the site where the wind is calculated.

Then, the Monin–Obukhov length scale $L$ is extracted from the mesoscale simulation and used to calculate the stability correction factor $C_U^{\mathrm{stability}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$ which is also a function of the point $\left(x\right)$ and wind direction ${\theta }_{\mathrm{ref}}^{\mathrm{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.

$$U_{\mathrm{target}}^{\mathrm{initial}}\left(x,t\right)=U_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)\times C_U^{\mathrm{stability}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$$

(4)

$${\sigma }_{\mathrm{target}}^{\mathrm{initial}}\left(x,t\right)={\sigma }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)$$

(5)

$${\theta }_{\mathrm{target}}^{\mathrm{initial}}\left(x,t\right)={\theta }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,t\right)$$

(6)

Finally, an advanced wake model is applied to calculate the velocity deficit $\Delta U\left(x,y,z,t\right)$, and added turbulence $\Delta \sigma \left(x,y,z,t\right)$ from wind turbines. The detail of the wake model is described in Section 2.3.

$$U_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)=U_{\mathrm{target}}^{\mathrm{initial}}\left(x,t\right)-\Delta U\left(x,t\right)$$

(7)

$$I_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)=\frac{\sqrt{{{\sigma }_{\mathrm{target}}^{\mathrm{initial}}}^2\left(x,t\right)+\Delta {\sigma }^2\left(x,t\right)}}{U_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)}$$

(8)

$${\theta }_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)={\theta }_{\mathrm{target}}^{\mathrm{initial}}\left(x,t\right)$$

(9)

where $U_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)$, $I_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)$ and ${\theta }_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)$ are the predicted mean wind speed, turbulence 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.

$$U_{\mathrm{target}}^{\mathrm{pred}}\left(x,t\right)=U_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)\times C_U^{\mathrm{terrain}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)\times C_U^{\mathrm{stability}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)-\Delta U\left(x,t\right)$$

(10)

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_U^{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ and $C_U^{\mathrm{stability}}$. 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.

$$\frac{\partial \rho {\tilde{u}}_j}{\partial x_j}=0$$

(11)

$$\frac{\partial \rho {\tilde{u}}_i}{\partial t}+\frac{\partial \rho {\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=-\frac{\partial \tilde{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\tilde{u}}_i}{\partial x_j}\right)+\frac{\partial {\tau }_{ij}}{\partial x_j}+f_{u,i}$$

(12)

where ${\tilde{u}}_i$ is the filtered wind speed in the x_i direction. $\rho$ is the air density, $\tilde{p}$ is the filtered pressure and $\mu$ is the molecular viscosity of the air.

${\tau }_{ij}$ is the SGS stresses which is modeled as follows:

$${\tau }_{ij}=2{\mu }_t{\tilde{S}}_{ij}+\frac{1}{3}{\tau }_{kk}{\delta }_{ij}$$

(13)

$${\tilde{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}\right)$$

(14)

where ${\mu }_t$ denotes the SGS turbulent viscosity, ${\tilde{S}}_{ij}$ is the rate of strain tensor for the resolved scale and ${\delta }_{ij}$ is the Kronecker delta.

Smagorinsky–Lilly is used for the SGS turbulent velocity:

$${\mu }_t=\rho L_s^2\left|\tilde{S}\right|=\rho L_s\sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}$$

(15)

$$L_s=\min \left(\kappa d, C_s{V}^{1/3}\right)$$

(16)

where $L_s$ denotes the mixing length for the sub-grid scales, $\kappa$ 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).

$$f_{u,i}=-\frac{F_{u,i}}{V_{grid}}=\frac{1}{2}\rho C_f\frac{1}{l_0}{\gamma }_0\sqrt{{\overline{u}}_j^2}{\overline{u}}_{i }, C_f=\frac{C_D}{{\left(1-{\gamma }_0\right)}^2}$$

(17)

where $c_f$ is a kind of drag coefficient representative, $l_0$ is the representative length scale of the obstacle and ${\gamma }_0$ is the packing density. The values of $c_f$, $l_0$ and ${\gamma }_0$ need to be set according to different types of obstacles. To consider the effect of vegetation, zero-displacement 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).

$${\tilde{C}}_U^{\mathrm{terrain}}\left(x,n_{\theta }\right)=\frac{U_{\mathrm{target}}^{\mathrm{terrain}}\left(x,n_{\theta }\right)}{U_{\mathrm{ref}}^{\mathrm{terrain}}\left(x_{\mathrm{ref},},n_{\theta }\right)}$$

(18)

$${\tilde{C}}_{\sigma }^{\mathrm{terrain}}\left(x,n_{\theta }\right)=\frac{{\sigma }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,n_{\theta }\right)}{{\sigma }_{\mathrm{ref}}^{\mathrm{terrain}}(x_{\mathrm{ref},}{n}_{\theta })}$$

(19)

$${\widetilde{\Delta \theta }}^{\mathrm{terrain}}\left(x,n_{\theta }\right)={\theta }_{\mathrm{target}}^{\mathrm{terrain}}\left(x,n_{\theta }\right)-{\theta }_{\mathrm{ref}}^{\mathrm{terrain}}(x_{\mathrm{ref}},n_{\theta })$$

(20)

in which ${\sigma }^{\mathrm{terrain}}$ is the fluctuating wind speed, $U^{\mathrm{terrain}}$ is the mean wind speed and ${\theta }^{\mathrm{terrain}}$ is the wind direction obtained from the CFD results and $n_{\theta }$ is the index number of direction sector. The correction factors for any direction, ${\theta }_{\mathrm{ref}}^{\mathrm{obs}}$, are obtained using a linear interpolation.

First, the direction sector of ${\theta }_{\mathrm{ref}}^{\mathrm{obs}}$ is obtained as shown in Equation (21).

$$n_{\theta }=\{\begin{array}{c}\lfloor \frac{{\theta }_{\mathrm{ref}}^{\mathrm{obs}}+11.25}{22.5}\rfloor {\theta }_{\mathrm{ref}}^{\mathrm{obs}}<348.75 \\ 0 \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}\end{array}$$

(21)

The linear interpolation is then performed as follows.

$$C_{\mathrm{U}}^{\mathrm{terrain}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)={\tilde{C}}_U^{\mathrm{terrain}}\left(x,n_{\theta }\right)+a\left({\tilde{C}}_U^{\mathrm{terrain}}\left(x,n_{\theta }+1\right)-{\tilde{C}}_U^{\mathrm{terrain}}\left(x,n_{\theta }\right)\right)$$

(22)

$$C_{\sigma }^{\mathrm{terrain}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)={\tilde{C}}_{\sigma }^{\mathrm{terrain}}\left(x,n_{\theta }\right)+a\left({\tilde{C}}_{\sigma }^{\mathrm{terrain}}\left(x,n_{\theta }+1\right)-{\tilde{C}}_{\sigma }^{\mathrm{terrain}}\left(x,n_{\theta }\right)\right)$$

(23)

$$\Delta {\theta }^{\mathrm{terrain}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)={\widetilde{\Delta \theta }}^{\mathrm{terrain}}\left(x,n_{\theta }\right)+a\left({\widetilde{\Delta \theta }}^{\mathrm{terrain}}\left(x,n_{\theta }+1\right)-{\widetilde{\Delta \theta }}^{\mathrm{terrain}}\left(x,n_{\theta }\right)\right)$$

(24)

where the proportional coefficient $a$ is defined in Equation (24).

$$a=\frac{{\theta }_{\mathrm{ref}}^{\mathrm{obs}}-{\theta }_{\mathrm{ref}}^{\mathrm{terrain}}\left(n_{\theta }\right)}{{\theta }_{\mathrm{ref}}^{\mathrm{terrain}}\left(n_{\theta }+1\right)-{\theta }_{\mathrm{ref}}^{\mathrm{terrain}}\left(n_{\theta }\right)}$$

(25)

2.2. 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. List of parameters in WRF.

Types	Parameter	Definition
Prognostic variables	$v\left(u,v,w\right)$	Covariant velocities
	$\varphi =gz$	Geopotential
	${\theta }_m$	Moist potential temperature
	$q_v$, $q_c$, $q_r$, $q_i$	mixing ratios for water vapor, cloud, rain, ice
Diagnostic variables	p	Pressure
	$P_{dh}$	Hydrostatic component of pressure for dry air
	$P_{ht}$	$P_h$ along the top boundary for dry air
	${\alpha }_d=1/\rho$	Reciprocal of dry atmospheric density
	$\theta$	Potential temperature
	${\mu }_d$	mass of dry air in the column
Constants	$g$	Gravitational acceleration
	$\gamma =C_p/C_v$	Specific heat ratio
	$p_0$	Reference pressure ${10}^5 \mathrm{Pa}$
	$R_d$	Gas constant for dry air
	$R_v$	Gas constant for water vapor
Predefined variables	$\eta$	Vertical coordinate
Forcing terms	$F_U$, $F_V F_W$, $F_{\mathrm{\Theta }}$, $F_{Q_m}$	External forces

.

The ARW equations are formulated using a terrain-following hydrostatic pressure vertical coordinate denoted by $\eta$.

$$\eta =\frac{P_{dh}-P_{dht}}{{\mu }_d}$$

(26)

where ${\mu }_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 = ${\mu }_{d}v$ and $\mathrm{\Theta }={\mu }_d\theta$. The flux equations are:

$$\frac{\partial U}{\partial t}+\left(\mathbf{\nabla }·Vu\right)+{\mu }_d\alpha \frac{\partial p}{\partial x}+\left(\frac{\alpha }{{\alpha }_d}\right)\frac{\partial p}{\partial \eta }\frac{\partial \varphi }{\partial x}=F_U$$

(27)

$$\frac{\partial V}{\partial t}+\left(\mathbf{\nabla }·Vv\right)+{\mu }_d\alpha \frac{\partial p}{\partial y}+\left(\frac{\alpha }{{\alpha }_d}\right)\frac{\partial p}{\partial \eta }\frac{\partial \varphi }{\partial y}=F_V$$

(28)

$$\frac{\partial W}{\partial t}+\left(\mathbf{\nabla }·Vw\right)+g\left[\left(\frac{\alpha }{{\alpha }_d}\right)\frac{\partial p}{\partial \eta }-{\mu }_d\right]=F_W$$

(29)

The governing thermodynamics equation, conservation of mass, geopotential time derivative and the equation of state are, respectively, shown in Equations (30)–(33).

$$\frac{\partial \mathrm{\Theta }}{\partial t}\mathbf{+}\nabla ·V\theta =F_{\mathrm{\Theta }}$$

(30)

$$\frac{\partial {\mu }_d}{\partial t}+\mathbf{\nabla }·V=0$$

(31)

$$\frac{\partial \Phi }{\partial t}+\frac{1}{{\mu }_d}\left[\left(V·\mathbf{\nabla }\varphi \right)-gW\right]=0$$

(32)

$$p=p_0{\left(\frac{R_d{\theta }_m}{p_0{\alpha }_d}\right)}^{\gamma }$$

(33)

$$\frac{\partial Q_m}{\partial t}+\left(\mathbf{\nabla }·\mathbf{V}{q}_m\right)=F_{Q_m}$$

(34)

where $\alpha$ is the inverse density taking the whole parcel into account $\alpha ={\alpha }_d{\left(1+q_v+q_c+q_r+q_i\right)}^{-1}$ and $Q_m={\mu }_d{q}_m; q_m=q_v, q_c, q_i.$

Additionally, ${\theta }_m=\theta \left(1+\left(\frac{R_v}{R_d}\right)q_v\right)\approx \theta \left(1+1.61q_v\right)$.

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].

$$\frac{\partial \theta }{\partial t}=F\left(\theta \right)+G_{\theta }\left(\widehat{{\theta }_0}-\theta \right)$$

(35)

where $\theta$ is the prognostic variable, $F\left(\theta \right)$ represents the normal tendency terms due to physics, the nudging parameter, $G_{\theta }$ is the timescale controlling the nudging strength and $\widehat{{\theta }_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].

$$u\left(z\right)=\frac{u_{*}}{\kappa }\left[ \ln \left(\frac{z}{z_0}\right)-\psi \left(\frac{z}{L}\right)\right]$$

(36)

where $u_{*}$ is the friction velocity, $\kappa$ 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 $\psi \left(\frac{z}{L}\right)$ is defined in Equation (37).

$$\psi \left(\frac{z}{L}\right)=\{\begin{array}{cc}-15.14& \left(7\le \frac{z}{L}\right)\\ -0.7\left(\frac{z}{L}\right)-0.75\left(\frac{z}{L}-\frac{5}{0.35}\right)e^{-0.35\frac{z}{L}}-\frac{0.75\times 5}{0.35}& \left(0.5<\frac{z}{L}<7\right)\\ -5\left(\frac{z}{L}\right)& \left(0<\frac{z}{L}\le 0.5\right)\\ 0& \left(\frac{z}{L}=0\right)\\ \begin{array}{c}2\ln \left(\frac{1+x}{2}\right)+\ln \left(\frac{1+x^2}{2}\right)-2{\tan }^{-1}\left(x\right)+\frac{\pi }{2}\\ x={\left[1-\left(16\frac{z}{L}\right)\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4 $}\right.}\end{array}& \left(\frac{z}{L}<0\right) \end{array}$$

(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:

$$C_{U\_t}^{\mathrm{stability}}\left(x,t\right)=\frac{\frac{\left[\ln \left(\frac{z}{z_0}\right)-\psi \left(\frac{z}{L\left(x,y,t\right)}\right)\right]}{\left[\ln \left(\frac{z_{\mathrm{ref}}}{z_0}\right)-\psi \left(\frac{z_{\mathrm{ref}}}{L\left(x,y,t\right)}\right)\right]}}{\frac{\left[\ln \left(\frac{z}{z_0}\right)\right]}{\left[\ln \left(\frac{z_{\mathrm{ref}}}{z_0}\right)\right]}}$$

(38)

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_{U\_t}^{\mathrm{stability}}\left(x,t\right)$ is defined assuming the terrain is flat. The value $C_{U\_t}^{\mathrm{stability}}\left(x,t\right)$ becomes unity for the same height ($z=z_{\mathrm{ref}}$) or neutral ($\psi =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).

$${\tilde{C}}_U^{\mathrm{stability}}\left(x,n_{\theta }\right)=\frac{\frac{\left[\ln \left(\frac{z}{z_0}\right)-{\psi }_{\mathrm{eq}}\left(x,n_{\theta }\right)\right]}{\left[\ln \left(\frac{z_{\mathrm{ref}}}{z_0}\right)-{\psi }_{\mathrm{eq}}\left(x,y,z_{\mathrm{ref}},n_{\theta }\right)\right]}}{\frac{\left[\ln \left(\frac{z}{z_0}\right)\right]}{\left[\ln \left(\frac{z_{\mathrm{ref}}}{z_0}\right)\right]}}$$

(39)

where ${\psi }_{\mathrm{eq}}\left(x,n_{\theta }\right)$ can be calculated as shown in Equation (39).

$${\psi }_{\mathrm{eq}}\left(x,n_{\theta }\right)=\frac{{{\displaystyle \sum }}^{ }{U}_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)\psi \left(\frac{z}{L\left(x,y,t\right)}\right)}{{{\displaystyle \sum }}^{ }{U}_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)}$$

(40)

$U_{\mathrm{ref}}^{\mathrm{obs}}\left(t\right)$ is the observed wind speed in that direction sector at time t. Once ${\psi }_{\mathrm{eq}}$ is calculated for each direction sector, the equivalent value of 1/$L_{\mathrm{eq}}$ for this wind direction can be calculated inversely which can explain the effect of stability for this wind direction and expressed as:

$$\frac{1}{L_{\mathrm{eq}}\left(x,y,n_{\theta }\right)}=\frac{{\psi }^{-1}\left({\psi }_{\mathrm{eq}}\left(x,n_{\theta }\right)\right)}{z}$$

(41)

To obtain $C_U^{\mathrm{stability}}\left(x,{\theta }_{\mathrm{ref}}^{\mathrm{obs}}\right)$ at any direction ${\theta }_{\mathrm{ref}}^{\mathrm{obs}}$, the ${\tilde{C}}_U^{\mathrm{stability}}\left(x,n_{\theta }\right)$ is linearly interpolated in the similar way as explained in Section 2.1.

2.3. 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. Summary of the Ishihara–Qian wake model.

Wake Model	Formula	Parameters
Wake width	$\frac{\sigma }{D}=k^{*}\frac{x}{D}+{\epsilon }^{*}, D_w=4\sqrt{2\ln 2}\sigma$	$k^{*}=0.11C{\prime }_T^{1.07}{I}_a^{0.20}$  ${\epsilon }^{*}=0.23{C^{\prime }}_T-0.25I_a^{0.17}$  ${C^{\prime }}_T=C_T\left(U_h,\gamma \right)\cos \mathsf{\gamma }, I_a>0.03$
Velocity Deficit	$\frac{\Delta U\left(x,y,z\right)}{U_h}=\frac{1}{{\left\{a+b·\frac{x}{D}+c{\left(1+\frac{x}{D}\right)}^{-2}\right\}}^2}\exp \left(-\frac{r^2}{2{\sigma }^2}\right)$  $r=\sqrt{{\left(y-y_d\right)}^2-{\left(z-H\right)}^2}$	$a=0.93C{\prime }_T^{-0.75}{I}_a^{0.17}$  $b=0.42C{\prime }_T^{0.6}{I}_a^{0.2}$  $c=0.15C{\prime }_T^{-0.25}{I}_a^{-0.7}$  $d=2.3C{\prime }_T^{-1.2}, e=1.0I_a^{0.1}$
Added Turbulence	$\frac{\Delta \sigma \left(x,y,z\right)}{U_h}=\frac{1}{d+e·\frac{x}{D}+f{\left(1+\frac{x}{D}\right)}^{-2}}\{k_1\exp \left(-\frac{{\left(r-\raisebox{1ex}{$D$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}^2}{2{\sigma }^2}\right)+k_2\exp \left(-\frac{{\left(r+\raisebox{1ex}{$D$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}^2}{2{\sigma }^2}\right)-\Delta I_a$  $\Delta I_a=\{\begin{array}{c}0 \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}\\ I_a{\sin }^2\left(\pi \frac{H-z}{D}\right){\cos }^2\left(\pi \frac{y}{D}\right) (0\le z<H, \left|y\right|\le D)\end{array}$	$f=0.7C_T^{\prime -3.2}{I}_a^{-0.45}$  $k_1\{\begin{array}{c}{\cos }^2\left(\frac{\pi }{2}\left(\frac{r}{D}-0.5\right)\right) \frac{r}{D}\le 0.5\\ 1 \frac{r}{D}>0.5\end{array}$  $k_1\{\begin{array}{c}{\cos }^2\left(\frac{\pi }{2}\left(\frac{r}{D}+0.5\right)\right) \frac{r}{D}\le 0.5\\ 0 \frac{r}{D}>0.5\end{array}$

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.

[16]. This analytical Gaussian-based model provides the three-dimensional wake characteristics including wake width, velocity deficit and added turbulence.

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 ${\sigma }_u^{ }$ over the rotor divided by the rotor onset wind speed $U_h$ [16].

$$U_h=\frac{1}{A}{{\displaystyle \int }}_{\mathrm{rotor}}^{ }U\left(x,y,z\right)dA$$

(42)

$$I_a=\frac{1}{A{U}_h}\sqrt{{{\displaystyle \int }}^{ }{\sigma }_u^2\left(x,y,z\right)dA}$$

(43)

Using these values, the effect of wake on the flow field, $\Delta U$ and $\Delta \sigma$, 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 wind farm 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 the assumption of homogenous flow in their wind vector reconstruction as shown in [27]. This error is estimated and compensated using the CFD simulation as explained in Appendix A. The position of the wind turbines and the location of the lidar and Met mast sites are shown in Figure 2b.

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. These wake-affected directions in which wind turbines exist within 10$D$, where $D=62 \mathrm{m}$ is the rotor diameter in this wind farm, are shaded in the figure.

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. CFD simulation domain and mesh setting.

Domain and Mesh Setting
Minimum horizontal mesh size	7.5 m
Maximum horizontal mesh size	240 m
Minimum mesh domain size	800 m (Diameter)
Horizontal mesh stretching ratio	1.15
Minimum vertical mesh size	5 m
Vertical mesh stretching ratio	1.1

. 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.

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. Description of mesoscale model and computational conditions.

	Domain 1	Domain 2	Domain 3
Simulation period	November 2018–June 2019
Center point	45° N 143° E
Horizontal resolution	18 km ($100\times 100$ grids)	6 km ($130\times 154$ grids)	2 km ($100\times 100$ grids)
Vertical resolution	45 levels (Surface to 50 hPa)
Time step	72 s	24 s	8 s
Spin-up	Two days
Boundary condition	Era 5 0.25° × 0.25° hourly [31]
Physics option	Ferrier (new Eta) microphysics scheme [32] Rapid radiative transfer model [33] Dudhia scheme [34] Unified Noah land surface model [35] Mellor–Yamada–Janjic (Eta) TKE level 2.5 scheme [36] Betts–Miller–Janic scheme (Except for domain 3) [37]

. 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.

In this study, $G_{\theta }$ 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) 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 the minimum requirement of the mean error of 5%.

3.3. 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. Classification of atmospheric stability according to Monin–Obukhov length scale.

Atmospheric Stability	Monin–Obukhov Length
Very stable	10 < L < 50
Stable	50 < L < 200
Neutral	|L| > 200
Unstable	−200 < L < −100
Very unstable	−100 < L < −50

. 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 equivalent value of 1/$L_{\mathrm{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] and 4 m/s [39] are filtered. Figure 12 shows the equivalent value of 1/$L_{\mathrm{eq}}$ by using these criteria. The harmonic average of L proposed by Duran et al. [14] is also plotted in Figure 12.

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 Figure 14 and Figure 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.

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.

• 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.
• A stability correction factor is proposed using the Monin–Obukhov length scale obtained from the mesoscale model to evaluate the effect of atmospheric stability on the mean wind speed. The proposed equivalent Monin–Obukhov length scale gives better agreement with the measurements than the conventional harmonic average of Monin–Obukhov length scale.
• The effects of topography, atmospheric stability, and turbine wake on the vertical profile of wind speed and turbulence intensity are also investigated. The consideration of all three effects gives better agreement with the measurements by the vertical lidar.