1. Introduction
Extended wind farms lead to the formation of a fully developed wind-turbine array boundary layer in which the power extraction from the wind farm originates from vertical energy transport from the flow aloft (Stevens & Meneveau Reference Stevens and Meneveau2017; Porté-Agel, Bastankhah & Shamsoddin Reference Porté-Agel, Bastankhah and Shamsoddin2020). As wind-turbine arrays are covering increasingly larger areas, the fully developed wind farm flow has been extensively studied using wind-tunnel experiments (Chamorro, Arndt & Sotiropoulos Reference Chamorro, Arndt and Sotiropoulos2011; Bossuyt, Meneveau & Meyers Reference Bossuyt, Meneveau and Meyers2017; Segalini & Chericoni Reference Segalini and Chericoni2021), large-eddy simulations (LES, Calaf, Meneveau & Meyers Reference Calaf, Meneveau and Meyers2010; Lu & Porté-Agel Reference Lu and Porté-Agel2011; Goit & Meyers Reference Goit and Meyers2015) and analytical models (Frandsen Reference Frandsen1992; Calaf et al. Reference Calaf, Meneveau and Meyers2010; Abkar & Porté-Agel Reference Abkar and Porté-Agel2013). These studies have deepened our fundamental understanding of the turbulent energy transport in extended wind farms (Meneveau Reference Meneveau2019), which is relevant to accurately representing the impact of large wind farms in climate model simulations (Keith et al. Reference Keith, DeCarolis, Denkenberger, Lenschow, Malyshev, Pacala and Rasch2004; Wang & Prinn Reference Wang and Prinn2010; Li et al. Reference Li, Kalnay, Motesharrei, Rivas, Kucharski, Kirk-Davidoff, Bach and Zeng2018) and models that capture regional meteorology (Baidya-Roy, Pacala & Walko Reference Baidya-Roy, Pacala and Walko2004). Note that two main types of wind farm parametrizations are used in weather and climate models (Porté-Agel et al. Reference Porté-Agel, Bastankhah and Shamsoddin2020): In the first type, the wind farm is considered as an increased surface roughness length; and in the second type, the wind farm is considered as an elevated sink of momentum. While the former was used early both in climate and mesoscale models, the latter is becoming more popular in mesoscale models.
The most popular analytical method to describe fully developed wind-turbine arrays is the top-down approach, which models the impact of a wind farm on the atmospheric flow above by determining the effective surface roughness height that corresponds to the wind farm. It has been inspired by the impact of plant canopies on atmospheric flows and has been considered to model wind farm dynamics since the pioneering studies of Templin (Reference Templin1974) and Newman (Reference Newman1977). In particular, Frandsen (Reference Frandsen1992) obtained an explicit formulation for effective surface roughness by assuming that the vertical profile of the horizontally averaged wind speed can be split into two logarithmic layers, one below and one above the turbine hub height. The model relates the corresponding friction velocities by balancing the momentum fluxes and assuming that the wind speed is continuous at hub height. The existence of the two logarithmic layers, with a wake layer in between, has been confirmed by the wind farm simulation of Calaf et al. (Reference Calaf, Meneveau and Meyers2010). They also introduced the concept of wake eddy viscosity to represent the wake layer and showed that it gives improved predictions compared with the model of Frandsen (Reference Frandsen1992). In the models of Frandsen (Reference Frandsen1992) and Calaf et al. (Reference Calaf, Meneveau and Meyers2010), the horizontally averaged wind speed at hub height is considered to be representative of the velocity at the turbine locations. However, the appropriateness of this assumption depends on the wind farm layout. To account for the local variations in the wind speed at hub height, Yang, Kang & Sotiropoulos (Reference Yang, Kang and Sotiropoulos2012) introduced the inhomogeneity factor 
$\beta$, which is the ratio of the velocity at the turbine location over the horizontally averaged velocity. Using one-dimensional mass and momentum conservation, they determined 
$\beta$ by approximating the mean streamwise velocity distributions at hub height. Stevens, Gayme & Meneveau (Reference Stevens, Gayme and Meneveau2015, Reference Stevens, Gayme and Meneveau2016b) introduced an ‘effective wake area coverage’ 
$w_f$ in their generalized coupled wake boundary-layer model, which is equivalent to setting 
$\beta ^{2} = {1/w_{f}}$. More recently, Zhang et al. (Reference Zhang, Ge, Liu and Yang2021) used a revised model of Jensen (Reference Jensen1983) to capture the effect of the turbine positioning on the velocity distribution in the wind farm and the model of Calaf et al. (Reference Calaf, Meneveau and Meyers2010) to simulate the effect of the vertical kinetic energy flux on the wind farm performance. These two models were coupled such that 
$\beta$ can be solved iteratively. Note that these models are all explicitly based on the logarithmic law. In contrast, Nishino & Dunstan (Reference Nishino and Dunstan2020) derived a two-scale momentum theory directly from the law of momentum conservation without using the logarithmic law, so that the theory may account for the effect of large-scale motions of the atmospheric boundary layer (ABL) in a time-dependent (rather than statistical) manner.
The top-down models are derived from truly neutral ABLs driven by a mean pressure gradient. However, ABLs are seldom truly neutral but often capped by a stably stratified free atmosphere, which is known as conventionally neutral ABLs (Zilitinkevich & Esau Reference Zilitinkevich and Esau2002). The mean pressure gradient is balanced by the Coriolis force associated with the geostrophic wind when the flow is in geostrophic equilibrium. This leads to the geostrophic drag law (GDL) for flat terrain, which was originally derived by Rossby & Montgomery (Reference Rossby and Montgomery1935) from a turbulence closure model and later by Blackadar & Tennekes (Reference Blackadar and Tennekes1968) from a more general asymptotic similarity. For extended wind farms, the equivalent wind farm surface roughness height is commonly used in conjuncture with the GDL to model the effect of the geostrophic wind forcing (Frandsen et al. Reference Frandsen, Barthelmie, Pryor, Rathmann, Larsen, Højstrup and Thøgersen2006; Calaf et al. Reference Calaf, Meneveau and Meyers2010; Meyers & Meneveau Reference Meyers and Meneveau2012),
\begin{equation} \left( \frac{\kappa G}{u_{*,2}} \right)^{2} = \left[ \ln \left( \frac{u_{*,2}}{|f| z_{0,2}} \right) - A \right]^{2} + B^{2}, \end{equation}
where 
$\kappa =0.4$ is the von Kármán constant, 
$G$ is the geostrophic wind speed, 
$u_{*,2}$ is the wind farm friction velocity, 
$z_{0,2}$ is the equivalent wind farm surface roughness, 
$f=2\varOmega \sin \phi$ is the Coriolis parameter, 
$\varOmega$ is the angular speed of the Earth and 
$\phi$ is the latitude. Typical values for the GDL coefficients 
$A$ and 
$B$ reported from meteorological observations are 
$A=1.3\unicode{x2013}1.9$ and 
$B=4.4\unicode{x2013}4.9$ (Hess & Garratt Reference Hess and Garratt2002). It has been shown that, for conventionally neutral ABLs, the coefficients 
$A$ and 
$B$ only depend on the Zilitinkevich number 
$Zi = N/f$ (Esau Reference Esau2004; Kadantsev, Mortikov & Zilitinkevich Reference Kadantsev, Mortikov and Zilitinkevich2021; Liu, Gadde & Stevens Reference Liu, Gadde and Stevens2021a,Reference Liu, Gadde and Stevensb; Liu & Stevens Reference Liu and Stevens2022). Here 
$N=\sqrt {g \varGamma /\theta _0}$ is the free-atmosphere Brunt–Väisälä frequency, 
$g$ is the gravity acceleration, 
$\varGamma$ is the free-atmosphere lapse rate and 
$\theta _0$ is the reference potential temperature. However, it is not a priori clear whether the similarity properties of the GDL apply to atmospheric flow above extended wind farms and, thus, whether a functional dependence that only depends on the Zilitinkevich number can describe the GDL coefficients 
$A$ and 
$B$.
To account for the effect of free-atmosphere stratification, Abkar & Porté-Agel (Reference Abkar and Porté-Agel2013) modified the model of Frandsen (Reference Frandsen1992) by adding a term 
$a_u N z$ to the logarithmic laws below and above hub height, where 
$a_u\approx 0.3$ is an empirical constant and 
$z$ is the altitude. The model system was closed by assuming the wind speed profile to be valid up to the boundary-layer height, where the wind speed equals the geostrophic wind speed. In agreement with their LES results, the model predicts that the power extracted by each turbine decreases with increasing stability of the free atmosphere. However, this model can not distinguish different wind farm layouts nor accurately capture the wind speed profile as the additional term 
$a_u N z$ will increase the wind speed at the hub height while the turbines impose an additional drag at this height. Also Allaerts & Meyers (Reference Allaerts and Meyers2015) proposed a simple analytical model to obtain the wind farm power output for the fully developed regime. With the assumption that the wind farm friction velocity is much larger than the surface friction velocity, the authors derived an approximate relation between the wind farm power output 
$P$ and the geostrophic drag coefficient 
$u_{*,2}/G$. The latter was determined by further assuming that the GDL of (1.1) was still valid above extended wind farms. Their model system is closed by estimating the equivalent wind farm surface roughness height 
$z_{0,2}$ from previous models (e.g. Frandsen Reference Frandsen1992; Calaf et al. Reference Calaf, Meneveau and Meyers2010) and determining explicit expressions for the GDL coefficients 
$A$ and 
$B$ with some specific eddy viscosity profiles (e.g. Csanady Reference Csanady1974; Nieuwstadt Reference Nieuwstadt1983). The resultant model is thus a function of the ABL height and the Rossby number. However, as remarked previously, the validity of the GDL above extended wind farms is not a priori clear and needs validation. In addition, as pointed out by the authors, the model of Allaerts & Meyers (Reference Allaerts and Meyers2015) is not valid when the wind farm thrust coefficient becomes very small, as in this case the wind farm friction velocity is no longer much larger than the surface friction velocity.
In the present work, we use the LES of wind farms to establish that insight from the GDL applies to fully developed wind-turbine array boundary layers in conventionally neutral ABLs. After providing parametrization of the GDL coefficients 
$A$ and 
$B$ for extended wind farms, we combine this updated GDL in conjuncture with an adjusted version of the model of Frandsen (Reference Frandsen1992) to model the wind farm production in conventionally neutral ABLs. Using this model, we investigate the effects of free-atmosphere lapse rate and turbine thrust coefficient on the wind farm power output and compare the results with LES data. The good agreement between the model prediction and LES data confirms the validity of the present model.
The organization of the paper is as follows. In § 2 the simulation method is explained. In § 3 the GDL for extended wind farms is parameterized based on the LES data. In § 4 we construct an analytical solution of the fully developed wind farm in conventionally neutral ABLs. In §§ 5 and 6 we investigate the effects of free-atmosphere lapse rate and turbine thrust coefficient, respectively. Finally, § 7 concludes the paper.
2. Large-eddy simulations of wind farms
We use an LES to simulate the conventionally neutral ABL flow over an infinite wind farm. In LES all turbulent structures larger than the filter scale are resolved, and the contributions of all smaller scales are parameterized using subgrid scale models. We integrate the spatially filtered governing equations (Albertson Reference Albertson1996; Liu & Stevens Reference Liu and Stevens2021a)
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} = 0, \end{gather}$$
$$\begin{gather}\partial_t \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{\omega}} \times \tilde{\boldsymbol{u}} = \boldsymbol{f}_{wt} + f \boldsymbol{e}_z \times (\boldsymbol{G} - \tilde{\boldsymbol{u}}) + g ( \tilde \theta / \langle \tilde \theta \rangle - 1) \boldsymbol{e}_z - \boldsymbol{\nabla} \tilde{p} / \rho - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$
$$\begin{gather}\partial_t \tilde{\theta} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\theta} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{q}}, \end{gather}$$
where 
$\tilde {{\cdot }}$ denotes the spatial filtering, 
$\langle \cdot \rangle$ represents the horizontal averaging, 
$\boldsymbol{e}_z$ is the unit vector in the vertical direction, 
$\tilde{\boldsymbol {u}}$ is the velocity, 
$\tilde {\boldsymbol {\omega } }= \boldsymbol {\nabla } \times \tilde {\boldsymbol {u}}$ is the vorticity, 
$\tilde p$ is the modified pressure departure from equilibrium, 
$\rho$ is the fluid density, 
$\boldsymbol {f}_{wt}$ is the force exerted by the wind turbines, 
$\tilde \theta$ is the potential temperature, 
$g$ is the gravity acceleration, 
$\boldsymbol {G}$ is the geostrophic wind velocity, 
$\boldsymbol {\tau }$ is the deviatoric part of the subgrid scale shear stress 
$\widetilde {\boldsymbol {u} \boldsymbol {u}} - \tilde {\boldsymbol {u}} \, \tilde {\boldsymbol {u}}$ and 
$\boldsymbol {q} = \widetilde {\boldsymbol {u} \theta } - \tilde {\boldsymbol {u}} \tilde \theta$ is the subgrid scale heat flux. Viscous terms in the governing equations are neglected as the Reynolds number is very high. The subgrid scale shear stress and heat flux are parameterized by the anisotropic minimum dissipation model (Abkar & Moin Reference Abkar and Moin2017; Gadde, Stieren & Stevens Reference Gadde, Stieren and Stevens2021).
The computational grid is uniform in the horizontal and vertical directions and staggered in the wall-normal direction. The first vertical velocity grid plane is located at the ground, while the first horizontal velocities grid plane is located at half a vertical grid distance above the ground. We use a pseudo-spectral discretization, and, thus, periodic boundary conditions, in the horizontal directions and a second-order finite difference method in the vertical direction. At the top boundary, we enforce a zero vertical velocity and zero shear stress. At the ground, we employ the classic wall model based on the logarithmic law of the wall (Moeng Reference Moeng1984; Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005). Time integration is performed using a second-order Adams–Bashforth method, and the projection method is used to ensure that the velocity field is divergence free. A proportional-integral-derivative controller is used such that the wind at hub height is always directed in the same way to ensure the wind farm layout is not changed due to the change of flow conditions in all simulations (Sescu & Meneveau Reference Sescu and Meneveau2014; Liu & Stevens Reference Liu and Stevens2021a). In addition, a Rayleigh damping layer is applied in the top 25 % of the domain to reduce the effects of wind farm induced gravity waves (Klemp & Lilly Reference Klemp and Lilly1978; Gadde & Stevens Reference Gadde and Stevens2021; Liu & Stevens Reference Liu and Stevens2021a).
The wind turbines are modelled using an actuator disk approach (Calaf et al. Reference Calaf, Meneveau and Meyers2010) in which the undisturbed upstream velocity 
$U_\infty$ is used to calculate the turbine force 
$F_{wt}$,
\begin{equation} F_{wt} ={-} \frac{1}{2}\rho C_T U_\infty^{2} \frac{\rm \pi}{4} D^{2}, \end{equation}
where 
$C_T$ is the thrust coefficient based on 
$U_\infty$ and 
$D$ is the turbine diameter. However, when a turbine operates in the wake of upstream turbines 
$U_\infty$ is not readily available. Therefore, 
$F_{wt}$ is actually calculated as (Calaf et al. Reference Calaf, Meneveau and Meyers2010)
\begin{equation} F_{wt} ={-} \frac{1}{2}\rho C_T' U_d^{2} \frac{\rm \pi}{4} D^{2}, \end{equation}
where 
$U_d$ is the disk-averaged velocity and 
$C_T'$ is the thrust coefficient based on 
$U_d$. From the momentum theory (e.g. Burton et al. Reference Burton, Sharpe, Jenkins and Bossanyi2001) it follows that 
$C_T=4a(1-a)$ and 
$C_T'=4a/(1-a)$, where 
$a$ is the axial induction factor. To avoid numerical instability, the turbine force is commonly filtered via a Gaussian smoothing kernel and this can result in a small overprediction of the wind-turbine production (Munters & Meyers Reference Munters and Meyers2017). It has been shown by Shapiro, Gayme & Meneveau (Reference Shapiro, Gayme and Meneveau2019) that this overprediction can be prevented by correcting the disk-averaged velocity 
$U_d$ in (2.5) by 
$U_d^{corr} = M U_d$, where 
$1/M=1+(C_T' \varDelta ) / (2 \sqrt {3{\rm \pi} } D)$ is the inverse correction factor and 
$\varDelta$ is the filter width. Thus, this correction method is always adopted in our simulations.
2.1. Conventionally neutral ABL simulations
The simulations are performed using an LES framework that has been validated extensively for flow in conventionally neutral ABLs (Liu et al. Reference Liu, Gadde and Stevens2021a,Reference Liu, Gadde and Stevensb) and wind farms (Stevens, Martínez-Tossas & Meneveau Reference Stevens, Martínez-Tossas and Meneveau2018; Liu & Stevens Reference Liu and Stevens2020, Reference Liu and Stevens2021b). Two sets of simulations are performed for aligned and staggered wind farm layouts. In the first set, we keep the thrust coefficient 
$C_T'=4/3$ constant to study the effect of the free-atmosphere lapse rate 
$\varGamma$ and latitude 
$\phi$ (table 1). In the second set, we keep the free-atmosphere lapse rate 
$\varGamma =4$ K km
$^{-1}$ constant to study the effect of the turbine thrust coefficient 
$C_T'$ and latitude 
$\phi$ (table 2). The computational domain size is 
$6\ {\rm km} \times 3\ {\rm km} \times 3\ {\rm km}$ in the streamwise, spanwise and vertical directions, and the corresponding grid resolution is 
$256\times 256\times 384$. Here streamwise and spanwise are relative to the mean wind direction at the hub height. We consider wind turbines with a rotor diameter 
$D=100$ m and hub height 
$z_h=100$ m. The turbine spacing normalized by the turbine diameter in the streamwise and spanwise directions is 
$s_x=s_y=6$, and, hence, the number of turbines is 
$50$ in each simulation. The initial potential temperature profile is 
$\theta (z) = \theta _0 + \varGamma z$ (Abkar & Porté-Agel Reference Abkar and Porté-Agel2013), where 
$\theta _0 = 300$ K is the potential temperature at the surface, 
$z$ is the altitude and 
$\varGamma = 0.2\sim 10$ K km
$^{-1}$. The initial velocity field is 
$\boldsymbol {u} = G \boldsymbol {e}_x$, where 
$G=12$ m s
$^{-1}$ is the geostrophic wind speed and 
$\boldsymbol {e}_x$ is the unit vector in the streamwise direction. To spin-up turbulence, small random perturbations are added to the initial fields of 
$\boldsymbol {u}$ and 
$\theta$ near the surface (
$z \le 100$ m). The surface roughness length is 
$z_{0,1}= 10^{-4}$ m, which is representative for offshore conditions, and the latitude is 
$\phi =30^{\circ }, 50^{\circ }$ and 
$80^{\circ }$. The statistics are averaged for the dimensionless time 
$ft\in [10, 10+2{\rm \pi} ]$, where the flow has reached its quasi-equilibrium state. In this state the mean potential temperature is nearly constant up to the capping inversion height (Abkar & Porté-Agel Reference Abkar and Porté-Agel2013). In the inversion layer there is a strong potential temperature gradient, which recovers to the free atmospheric lapse rate above in the free atmosphere. In total, 66 high-resolution simulations are performed and statistics are summarized in tables 1 and 2.
Summary of the first set of simulations with thrust coefficient 
$C_T'=4/3$.
Summary of the second set of simulations with lapse rate 
$\varGamma =4$ K km
$^{-1}$.
2.2. Boundary-layer height
Figure 1 shows the normalized boundary-layer height 
$u_{*,2}^{2}/(fh)^{2}$ as a function of the 
$Zi$ number, which reveals clearly that the boundary-layer height decreases as the stability of the free atmosphere increases (see also table 1). Here 
$u_{*,2}=(u_{*,1}^{2}+0.5 c_{ft}' U_d^{2})^{1/2}$ is the wind farm friction velocity, with 
$u_{*,1}$ being the surface friction velocity and 
$c_{ft}'={\rm \pi} C_T'/(4 s_x s_y)$ the wind farm thrust coefficient, and 
$h$ is the boundary-layer height, which is determined as the height where the total shear stress 
$\tau =0.05 u_{*,2}^{2}$ (Abkar & Porté-Agel Reference Abkar and Porté-Agel2013). The expression that describes well the LES data is 
$u_{*,2}^{2}/(fh)^{2}=Zi/C_N^{2}$, with 
$C_N=1.61$ being an empirical constant. This is very similar to the situation over flat terrain (Zilitinkevich, Esau & Baklanov Reference Zilitinkevich, Esau and Baklanov2007; Zilitinkevich et al. Reference Zilitinkevich, Tyuryakov, Troitskaya and Mareev2012), indicating that the similarity as defined in the GDL applies to fully developed wind-turbine array boundary layers. However, we note that the boundary-layer height above a wind farm is significantly higher than over flat terrain as the friction velocity 
$u_{*,2}$ in the former is much higher than in the latter. In addition, we note that the inversion height 
$z_i$, which is determined as the height where the total heat flux reaches its minimum value (Abkar & Moin Reference Abkar and Moin2017), can also be described well by 
$u_{*,2}^{2}/(fz_i )^{2}=Zi/C_N'^{2}$, where the empirical constant 
$C_N'=1.45< C_N=1.61$. Thus, the inversion height 
$z_i$ is generally shallower than the ABL height 
$h$, and their ratio is approximately 
$z_i/h \approx 0.9$.
The 
$Zi$-dependence of the inverse square of the dimensionless boundary-layer height 
$u_{*,2}^{2}/(fh)^{2}$. Filled symbols: LES data of table 1 for aligned (circles) and staggered (diamonds) wind farms; open symbols: the corresponding LES data of table 2; solid line: model prediction of 
$u_{*,2}^{2}/(fh)^{2}=Zi/C_N^{2}$ with 
$C_N=1.61$.
Figure 2 shows the velocity distribution at hub height for aligned and staggered wind farms normalized by the geostrophic wind speed. It is demonstrated that the wind velocity in this plane is non-uniform due to the wind-turbine wakes. High-velocity wind speed channels are formed between the turbines for the aligned wind farm configuration. At the same time, for the staggered layout, the wakes have a relatively long distance to recover. As a result, the velocity at the turbine location 
$U_d$ compared with the horizontally averaged velocity at hub height 
$U_h$ is relatively high for the staggered wind farms and relatively low for the aligned wind farms. Because the power output 
$P$ is proportional to 
$U_d^{3} = (1-a)^{3} \beta ^{3} U_h^{3}$ rather than 
$U_h^{3}$, it is important to model 
$\beta$ to account for the local velocity variations compared with the horizontal mean.
Temporally averaged streamwise velocity at hub height normalized by the geostrophic wind speed for (a,c,e) aligned and (b,d, f) staggered wind farms. The free-atmosphere lapse rate is 
$\varGamma =4$ K km
$^{-1}$ and the latitude is (a,b) 
$\phi =30^{\circ }$, (c,d) 
$\phi =50^{\circ }$ and (e, f) 
$\phi =80^{\circ }$. The short lines denote the location of turbines.
3. Geostrophic drag law for extended wind farms
Figure 3 shows the 
$Zi$-dependence of the GDL coefficients
\begin{equation} A = \ln{\left(\frac{u_{*,2}}{f z_{0,2}}\right)} - \frac{\kappa G}{u_{*,2}} \cos \alpha_0, \quad B = \frac{\kappa G}{u_{*,2}} \sin \alpha_0, \end{equation}
where 
$z_{0,2}$ is determined by 
$z_{0,2}=2z_h \exp (-\kappa U_{2h}/u_{*,2})$, with 
$U_{2h}$ being the horizontally averaged wind speed at 
$z=2z_h$ (Calaf et al. Reference Calaf, Meneveau and Meyers2010; Zhang et al. Reference Zhang, Ge, Liu and Yang2021), and 
$\alpha _0$ is the total wind angle change across the entire boundary layer (see table 1). It is demonstrated that 
$A$ and 
$B$ are well parametrized by a function that depends only on 
$Zi$, similar to what is observed in the classical GDL for flow over flat terrain (Zilitinkevich & Esau Reference Zilitinkevich and Esau2005; Liu et al. Reference Liu, Gadde and Stevens2021a). This confirms that the similarity arguments apply for flow over extended wind farm arrays, even though their physics are different. Because the friction velocity 
$u_{*,2}$ is significantly increased compared with the flat terrain case, the value of 
$B$ is lower, and the value of 
$A$ is higher for the wind farm case than for the traditional flat terrain case. Here we parameterize 
$A$ and 
$B$ for fully developed wind-turbine arrays in conventionally neutral ABLs using a least-square fitting to the data as follows:
\begin{equation} A=1.54+0.18 \ln {Zi}, \quad B=1.74+0.011 Zi. \end{equation}The corresponding results are shown in figure 3, which confirm that this model indeed captures the LES data well.
The 
$Zi$-dependence of the (a) 
$A$ and (b) 
$B$ coefficients in the GDL. Filled symbols: LES data of table 1 for aligned (circles) and staggered (diamonds) wind farms; open symbols: the corresponding LES data of table 2; solid line: the least-square fitting using (3.2a,b); dashed line: coefficients over flat terrain (Liu et al. Reference Liu, Gadde and Stevens2021a).
It is worth pointing out that a complete description of 
$A$ and 
$B$ in extended wind farms involves many parameters, such as the hub height 
$z_h$, rotor diameter 
$D$, turbine spacing 
$(s_x,s_y)$, surface roughness height 
$z_{0,1}$, geostrophic wind speed 
$G$, turbine thrust coefficient, wind farm layout, latitude and free-atmosphere lapse rate. For simplicity, we fixed 
$(z_h,D,s_x,s_y,z_{0,1},G)$ in our simulations. As shown in figure 3, the data collapse well for different wind farm layouts and turbine thrust coefficients, and only depend on the Zilitinkevich number. However, these results cannot guarantee the values of 
$A$ and 
$B$ are unchanged when 
$(z_h,D,s_x,s_y,z_{0,1},G)$ are varied, which should be kept in mind for future studies.
4. A new top-down model for conventionally neutral ABLs
Similar to the models of Frandsen (Reference Frandsen1992) and Abkar & Porté-Agel (Reference Abkar and Porté-Agel2013), we assume that the vertical profiles of the horizontally averaged wind speed 
$U$ is characterized by two logarithmic velocity regions, i.e. one below hub height and one above. Figure 4 shows a sketch of the corresponding wind speed profile, which is modelled as
\begin{equation} U =
\left\{\begin{array}{@{}ll} \displaystyle\dfrac{u_{*,
1}}{\kappa} \ln \left( \dfrac{z}{z_{0, 1}} \right) - a_u
\beta^{2} u_{*,1} \dfrac{z}{z_h}, & z \le z_h,\\
\displaystyle\dfrac{u_{*, 2}}{\kappa} \ln \left(
\dfrac{z}{z_{0, 2}} \right), & z \ge z_h.
\end{array}\right. \end{equation}
Here 
$a_u=a_u(C_T')$ is an empirical constant that can be determined from LES data by setting 
$z=z_h$ in the first formula of (4.1). For the cases considered here, we find that 
$a_u \approx 4.3$ for 
$C_T'=4/3$; see figures 5(a) and 8(a). The factor 
$\beta =U_d/[(1-a)U_h]$, which accounts for the flow inhomogeneity at hub height, is incorporated into the first line of (4.1). To determine 
$\beta$ analytically, Zhang et al. (Reference Zhang, Ge, Liu and Yang2021) introduced a new coupled model that combines a revised model of Jensen (Reference Jensen1983) and the top-down model of Calaf et al. (Reference Calaf, Meneveau and Meyers2010). Here we adopted the model of Zhang et al. (Reference Zhang, Ge, Liu and Yang2021) and added ghost turbines (Lissaman Reference Lissaman1979; Ge et al. Reference Ge, Yang, Zhang and Zuo2021) to further account for the ground effects. We find that 
$\beta =0.973$ for an aligned wind farm and 
$\beta =1.102$ for a staggered wind farm. Figure 5(b) shows that these values agree within 
$4\,\%$ of the values determined from LES. There may be various reasons for the observed small differences. First of all, the model of Zhang et al. (Reference Zhang, Ge, Liu and Yang2021) neglects the effects of the Coriolis force and thermal stratification. Furthermore, the performance of the model may depend on the parameters such as 
$s_x,s_y,z_{0,1}$, the wind farm layout, the employed wake superposition method and the boundary-layer height parametrization.
A sketch showing the different parameters used in the present model.
The 
$Zi$-dependence of (a) the empirical constant 
$a_u$ and (b) the factor 
$\beta$ for aligned (circles) and staggered (diamonds) wind farms. Solid line: (a) 
$a_u=4.3$ and (b) analytical predictions for aligned (
$\beta =0.973$) and staggered (
$\beta =1.102$) wind farms, see details in text.
The wind speed profile should be continuous at hub height, and from (4.1) we thus obtain
\begin{equation} U_h = \frac{u_{*, 1}}{\kappa} \ln \left( \frac{z_h}{z_{0, 1}} \right) - a_u \beta^{2} u_{*,1} =\frac{u_{*, 2}}{\kappa} \ln \left( \frac{z_h}{z_{0, 2}} \right). \end{equation}
Furthermore, the surface friction velocity 
$u_{*,1}$, the wind farm friction velocity 
$u_{*,2}$ and the hub height wind speed 
$U_h$ can be related by the horizontally averaged momentum balance above and below the turbines. When neglecting the Coriolis effect in the turbine wake region, the relation reads
\begin{equation} u_{*, 2}^{2} = u_{*, 1}^{2} + \frac{1}{2} c_{ft} \beta^{2} U_h^{2}, \quad c_{ft} = \frac{{\rm \pi} C_T}{4 s_x s_y}. \end{equation}
In addition, 
$u_{*,2}$ and 
$z_{0,2}$ are related by the GDL (i.e. (1.1)) with 
$A$ and 
$B$ being parameterized by (3.2a,b). This means that the system of equations ((1.1), (4.2) and (4.3a,b)) can be solved for the three unknowns 
$(u_{*, 1}, u_{*, 2}, z_{0, 2})$.
Figure 6 compared the vertical profiles of the horizontally averaged velocity obtained from LES with the above corresponding profiles obtained from the model (i.e. (4.1)). Very close to the wall, the LES data agrees well with the logarithmic law that originates from the wall. However, around hub height the influence of the wake effects, which are modelled by the additional term in (4.1), becomes pronounced. The present model also predicts the wind speed profile above hub height very well. The only exception is the region around the inversion layer, where the effects of the Coriolis force and thermal stratification result in the formation of a low-level jet.
The normalized horizontally averaged wind speed 
$U/G$ profile for aligned (a,c,e) and staggered (b,d, f) wind farms. The free-atmosphere lapse rate 
$\varGamma =4$ K km
$^{-1}$ and the latitude (a,b) 
$\phi =30^{\circ }$, (c,d) 
$\phi =50^{\circ }$ and (e, f) 
$\phi =80^{\circ }$. Solid line: prediction of the present wind farm model (4.1); dashed-dotted line: the logarithmic law for the flat terrain very close to the ground. The vertical dashed lines denote the location of turbine rotor area.
5. Effect of free-atmosphere lapse rate
Figure 7(a,b) compares the 
$\varGamma$-dependence of the normalized hub height wind speed 
$U_h/G$ obtained from LES with the model prediction. The LES data shows that the normalized wind speed at hub height increases as the free-atmosphere lapse rate 
$\varGamma$ decreases or the latitude 
$\phi$ increases (see also table 1). In addition, it is revealed that the wind farm layout has a significant effect as 
$U_h$ for the aligned layout is about 
$6\sim 8\,\%$ higher than for the staggered layout. The reason is that staggered wind farms extract more energy from the flow and this increased drag results in a lower normalized velocity at hub height. Panels (c,d) show the corresponding normalized total wind farm power output per surface area 
$10^{3} P/G^{3}$, where 
$P= 0.5 c_{ft}' U_d^{3}$ is the wind farm power output per surface area, 
$c_{ft}'={\rm \pi} C_T'/(4 s_x s_y)$ is the normalized wind farm thrust coefficient and 
$U_d=(1-a)\beta U_h$ is the velocity at the turbine location. These panels show that the predicted power output agrees well with the LES data. As expected, the power output is slightly higher in the staggered wind farms than in the aligned ones.
The 
$\varGamma$-dependence of the (a,b) normalized wind speed at hub height 
$U_h/G$ and the corresponding (c,d) normalized wind farm power output per surface area for (a,c) aligned and (b,d) staggered wind farm. Lines: present model; filled symbols: LES data.
It should be pointed out that the differences in wind farm performance between aligned and staggered layouts are strongly dependent on the spanwise spacing (among other parameters). Thus, in the fully developed regime, a staggered layout does not always lead to a higher power production. For example, the present study focuses on a uniform turbine spacing of 
$6D$ in the streamwise and spanwise directions and finds a higher power output for a staggered layout than for an aligned layout. However, previous LES and wind-tunnel experiments (Stevens, Gayme & Meneveau Reference Stevens, Gayme and Meneveau2014, Reference Stevens, Gayme and Meneveau2016a; Wu & Porté-Agel Reference Wu and Porté-Agel2017; Bossuyt, Meneveau & Meyers Reference Bossuyt, Meneveau and Meyers2018) have shown that, for certain conditions with smaller spanwise spacings (
${\le }5D$), staggered and aligned layouts can lead to the same asymptotic power output in the fully developed regime.
A closer inspection of figure 7 reveals that the predictions for the staggered layout are generally better than the aligned layout. This is reasonable as the flow in the staggered wind farm is closer to the horizontally homogeneous condition than the aligned layout (Stevens et al. Reference Stevens, Gayme and Meneveau2016a), see figure 2, which is one of the underlying assumptions of the model. Furthermore, it should be appreciated that wind farm designs of extended wind farms would be closer to a staggered layout as it optimizes wind farm performance. The figure also shows that the power decreases as the free atmospheric lapse rate increases. The reason is that kinetic energy entrainment, which is the only source of kinetic energy that balances the power extracted by the turbines, decreases as the stability in the free atmosphere increases. As a result, less power can be extracted from the flow by the wind turbines (Abkar & Porté-Agel Reference Abkar and Porté-Agel2013, Reference Abkar and Porté-Agel2014). In addition, one notes that the predictions for the hub height velocity tend to be closer to the LES values than the corresponding power production estimates. The reason is that the power production scales as 
$\beta ^{3} U_h^{3}$, making the power more difficult to predict than just the velocity 
$U_h$.
6. Effect of turbine thrust coefficient
To examine the validity of the proposed theoretical framework, we consider the second set of simulations in which the free-atmosphere lapse rate 
$\varGamma =4$ K km
$^{-1}$ is fixed, while the latitude 
$\phi$, the wind farm layout and the thrust coefficient 
$C_T'$ are varied. In particular, the thrust coefficient 
$C_T'$ changes between 
$C_T' = 0.1$ and 
$C_T' = 2$, i.e. the Betz limit; see table 2. Figure 8(a) shows that the empirical constant 
$a_u$ in (4.1) increases from zero at 
$C_T'=0$ and asymptotically approaches 
$a_u \approx 4.3$ for 
$C_T'=4/3$. The lower limit is consistent with the physical assumptions as it indicates that the wake correction term in (4.1) disappears when wind turbines are turned off. The fact that 
$a_u$ asymptotes to a constant value indicates that the deviation from the logarithmic law below the turbine hub height caused by wake effects saturates for a sufficiently large thrust coefficient, say 
$C_T' \ge 4/3$, as the energy entrainment from above is limited. For the analytical predictions, the value of 
$a_u$ is parametrized as 
$a_u=4.3\tanh (2C_T')$, which is in good agreement with the LES data as shown in figure 8(a). Moreover, figure 8(b) shows the comparison of the analytical predictions and LES measurements of the flow inhomogeneity factor 
$\beta$ for different thrust coefficients, of which the overall agreement is good for both the aligned and staggered wind farms. It is also presented that the value of 
$\beta$ is nearly independent of the latitude, indicating that our assumption of neglecting the Coriolis force in the wake region is reasonable.
The 
$C_T'$-dependence of (a) the empirical constant 
$a_u$ and (b) 
$\beta$ for aligned (circles) and staggered (diamonds) wind farms. The solid line in (a) 
$a_u=4.3\tanh (2C_T')$ and (b) 
$\beta$ is obtained from model calculations; see § 4 for details.
Figure 9(a,b) shows that the normalized wind speed at hub height 
$U_h/G$ reduces monotonically as the thrust coefficient 
$C_T'$ increases. As staggered wind farms extract more kinetic energy from the ABL flow than aligned wind farms under the same conditions, 
$U_h/G$ is smaller for the staggered wind farms. Furthermore, 
$U_h/G$ is larger at higher latitudes. This phenomenon originates from the increasing Coriolis force in the free atmosphere that drives the conventionally neutral ABL flow. With increasing latitude 
$\phi$, the Coriolis parameter 
$f=2\varOmega \sin \phi$ increases, allowing for an increased normalized wind speed at hub height 
$U_h/G$. It is worthwhile to point out that increasing 
$C_T'$ also increases the drag coefficient 
$u_{*,2}/G$ and the equivalent wind farm surface roughness 
$z_{0,2}$ perceived by the mean flow aloft, a trend supported well by the results presented in table 2.
The 
$C_T'$-dependence of the (a,b) normalized wind speed at hub height 
$U_h/G$ and the corresponding (c,d) normalized power output per surface area 
$10^{3}P/G^{3}$ for (a,c) aligned and (b,d) staggered wind farms. Lines: present model; filled symbols: LES data.
Figure 9(c,d) shows the corresponding normalized wind farm power outputs per surface area. When considering a single turbine in idealized conditions, the optimal operating condition of the turbine corresponds to 
$C_T'=2$, corresponding to the Betz limit (Burton et al. Reference Burton, Sharpe, Jenkins and Bossanyi2001). However, this figure shows that the normalized production of fully developed wind farms reaches its maximum value at 
$C_T'\approx 4/3$. A similar phenomenon has also been observed for the fully developed wind farm in truly neutral ABLs by Goit & Meyers (Reference Goit and Meyers2015). The reason is that the power production 
$P/G^{3}$ is proportional to the thrust coefficient 
$C_T'$ and the cubed hub height velocity 
$(U_h/G)^{3}$. As the thrust coefficient 
$C_T'$ increases, the hub height velocity 
$U_h/G$ decreases due to the wake effects (see table 2). For a low thrust coefficient, the benefit of increasing 
$C_T'$ is more significant than the decrease of 
$U_h/G$, leading to increased power output. However, the benefit of increasing 
$C_T'$ is overwhelmed by the reducing hub height velocity for a high thrust coefficient.
7. Conclusions
Using LES of fully developed wind farms in conventionally neutral ABLs, we demonstrate that the normalized boundary-layer height and the GDL coefficients 
$A$ and 
$B$ can be parameterized as a function of just the Zilitinkevich number 
$Zi$. This confirms that the similarity arguments from the GDL apply for flow over extended wind farms, even though their physics are different. The coefficients 
$A$ and 
$B$ are changed with respect to flow over flat terrain as the wind farm friction velocity is higher than over flat terrain. The canopy and wind farm flow similarity suggests that these findings may be relevant for canopy flow.
Inspired by the classical work of Frandsen (Reference Frandsen1992), we develop a new top-down model that predicts the performance of fully developed wind farm arrays as a function of the geostrophic wind forcing in conventionally neutral ABLs. The model captures the effect of the free-atmosphere lapse rate 
$\varGamma$, the latitude 
$\phi$, the wind farm layout and the thrust coefficient 
$C_T$ on the wind farm performance. We demonstrate that the highest wind farm power output is obtained for the thrust coefficient well below the Betz limit. This means that the optimal wind farm performance is obtained when the wind farm is operated as a complete system and not when each turbine optimizes its performance. Furthermore, the model predicts that the normalized power production per surface area decreases with the free-atmosphere lapse rate, which is consistent with the finding of Abkar & Porté-Agel (Reference Abkar and Porté-Agel2013). This phenomenon originates from the decreasing kinetic energy entrainment with increasing stability in the free atmosphere. In contrast, the power production per surface area normalized by the cubed geostrophic wind speed increases with the latitude, which is consistent with the finding of Antonini & Caldeira (Reference Antonini and Caldeira2021) and results from the increasing Coriolis force in the free atmosphere with increasing latitude.
Future research will focus on developing a similar model to include the effect of thermal stratification in the boundary layer. To do this, two main adjustments may be needed. The first one is to include a thermal stratification correction to the initially assumed logarithmic profile (Sescu & Meneveau Reference Sescu and Meneveau2015). The second one is to extend the GDL approach towards thermally stratified flows (Tong & Ding Reference Tong and Ding2020; Kadantsev et al. Reference Kadantsev, Mortikov and Zilitinkevich2021).
Acknowledgements
We thank Profs. Nansheng Liu and Mingwei Ge for their insightful discussions. We acknowledge Dr. Srinidhi N. Gadde for plotting the image in the Graphical Abstract.
Funding
This work was supported by the Hundred Talents Program of the Chinese Academy of Sciences (CAS), the National Natural Science Foundation of China Grant (No. 11621202), the Anhui NARI Jiyuan Electric Power Grid Technology Co. Ltd. through the Joint Laboratory of USTC-NARI, the Shell-NWO/FOM-initiative Computational sciences for energy research of Shell and Chemical Sciences, Earth and Live Sciences, Physical Sciences, FOM and STW, and an STW VIDI Grant (No. 14868). NWO Domain Science sponsored this work to use the national computer facilities. We acknowledge PRACE for awarding us access to MareNostrum4 in Barcelona, Spain. This research was also supported by the advanced computing resources provided by the Supercomputing Center of the USTC.
Declaration of interests
The authors report no conflict of interest.
Open access
