J. Fluid Mech. (2024), vol. 985, A43, doi:10.1017/jfm.2024.282

A fast-running physics-based wake model for asemi-infinite wind farm
Majid Bastankhah1, †, Mohammad Mehdi Mohammadi2 , Charlie Lees1 ,Gonzalo Pablo Navarro Diaz2 , Oliver R.H. Buxton3 and Stefan Ivanell2
1 Department of Engineering, Durham University, Durham DH1 3LE, UK
2 Department of Earth Sciences, Uppsala University, Visby 62167, Sweden
3 Department of Aeronautics, Imperial College London, London SW7 2AZ, UK

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press(Received 15 September 2023; revised 26 January 2024; accepted 14 March 2024)

This paper presents a new generation of fast-running physics-based models to predict thewake of a semi-infinite wind farm, extending infinitely in the lateral direction but withfinite size in the streamwise direction. The assumption of a semi-infinite wind farm enables
concurrent solving of the laterally averaged momentum equations in both streamwise andspanwise directions. The developed model captures important physical phenomena suchas vertical top-down transport of energy into the farm, variable wake recovery rate due tothe farm-generated turbulence and also wake deflection due to turbine yaw misalignment
and Coriolis force. Of special note is the model’s capability to predict and shed lighton the counteracting effect of Coriolis force causing wake deflections in both positive
and negative directions. Moreover, the impact of wind farm layout configuration on theflow distribution is modelled through a parameter called the local deficit coefficient.
Model predictions were validated against large-eddy simulations extending up to 45 km
downstream of wind farms. Detailed analyses were performed to study the impacts ofvarious factors such as incoming turbulence, wind farm size, inter-turbine spacing andwind farm layout on the farm wake.
Key words: wakes

1. Introduction

Offshore wind is projected to experience rapid expansion in the coming decades, emerging
as a significant global renewable energy source (Veers et al. 2019). To achieve this goal,
many new offshore wind farms are anticipated to be erected in specific and promising† Email address for correspondence: majid.bastankhah@durham.ac.uk© The Author(s), 2024. Published by Cambridge University Press. This is an Open Access article,
distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/
licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the originalarticle is properly cited.

985 A43-1

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

M. Bastankhah and others
geographical areas, particularly in regions like the North Sea, where strong and consistentwinds are present. Consequently, the interaction among neighbouring offshore wind farms,
as their wakes affect each other, has become an essential and pressing subject of research.
Recent satellite images and field measurements have revealed that wakes of wind farms canlast for many kilometres (Christiansen & Hasager 2005; Nygaard & Christian Newcombe
2018; Ahsbahs et al. 2020). Significant power degradation and fatigue loads can thus
occur for a wind farm subject to wakes of adjacent wind farms (Stevens & Meneveau2017). Beyond technical complexities, interactions between adjacent wind farms may lead
to legal and financial disputes between operators of neighbouring facilities. As a result,accurate and reliable modelling of wind farm wake effects becomes of great importancefor optimising future wind farms in increasingly competitive offshore environments.
High-fidelity numerical simulations such as large-eddy simulation (LES) are powerfultools for modelling complex turbulent wake flows, offering detailed insights into flow
dynamics and wake interactions (Porté-Agel, Bastankhah & Shamsoddin (2020), andreferences therein). However, simulating a cluster of wind farms in congested areas suchas the North Sea with LES is computationally intensive and time-consuming, making it
impractical for real-time or large-scale studies. To address this challenge and enable more
efficient simulations, there is a clear demand for fast-running engineering wake models
striking a balance between accuracy and computational cost. Major advantages of thesemodels are their ease of use and low computational costs, allowing for quicker assessmentsof various scenarios and aiding in optimisation of wind farm layouts and real-time control.
Below, we attempt to classify the engineering wake models developed in the literature.
The typical method for modelling airflow distribution within wind farms involves
predicting the wake generated by each individual turbine. A superposition method is thenapplied to consider the combined impact of these wake effects. These individual wakemodels range mainly from top-hat models (Jensen 1983; Katić, Højstrup & Jensen 1986)
to Gaussian-type models (Bastankhah & Porté-Agel 2014). The Jensen top-hat model (alsoknown as the Park model) has been extended in recent works to account for variablewake recovery rate due to turbine-generated turbulence (Nygaard et al. 2020). Over time,
Gaussian wake models have also been refined and extended in several studies to moreaccurately describe the near-wake region (e.g. Keane et al. 2016; Shapiro et al. 2019;Blondel & Cathelain 2020; Schreiber, Balbaa & Bottasso 2020), to better capture wake
expansion and its asymmetric shape (e.g. Abkar & Porté-Agel 2015; Xie & Archer 2015;
Pedersen et al. 2022; Vahidi & Porté-Agel 2022) or to capture effects of yaw angle (e.g.
Bastankhah & Porté-Agel 2016; King et al. 2020; Bastankhah et al. 2022; Bay et al. 2023)
and wind veer (Abkar, Sørensen & Porté-Agel 2018; Mohammadi et al. 2022; Narasimhan,
Gayme & Meneveau 2022). Moreover, a variety of wake superposition methods exist,
aiming to model cumulative wake effects in wind farms (e.g. Lissaman 1979; Voutsinas,
Rados & Zervos 1990; Niayifar & Porté-Agel 2016; Zong & Porté-Agel 2020; Bastankhah
et al. 2021; Lanzilao & Meyers 2022). Some of these methods are solely empirical innature, while others have a foundation in flow physics. See Bastankhah et al. (2021)
for a detailed discussion of different wake superposition methods. This simple approachhas proven to be very useful in providing detailed information on the flow field withinsmall-sized wind farms and has been extensively used in wind farm layout optimisation
and real-time flow control; see the review of Meyers et al. (2022) and references therein.
However, this modelling approach cannot properly describe the interaction of wind farms
with the atmospheric boundary layer (ABL) which involves scales that are comparable tothe size of the entire wind farm or the ABL thickness. Most notably, these models fall shortof capturing the crucial vertical transport of kinetic energy from higher-altitude layers of
985 A43-2

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Fast-running physics-based wake model for a wind farm
the atmosphere into the wind farm/wind farm wake (Stevens & Meneveau 2017). This
becomes especially problematic as the size of wind farms grows, or if we seek information
about the wake of an entire wind farm several kilometres downstream.
Capturing large-scale wind farm physics may be more readily achieved using infinite
wind farm models. In this approach, the wind farm is assumed to be infinitely large in bothlateral and streamwise directions, and the whole wind farm is modelled as an area with anincreased aerodynamic surface roughness. Unlike single-wake modelling, this approachis able to capture the vertical transport of energy caused by turbulent fluxes, which is inbalance with the energy extracted by wind turbines in infinite wind farms. The interested
reader is referred to the seminal works of Frandsen (1992) and Calaf, Parlange & Meneveau(2011) and other subsequent studies (e.g. Frandsen et al. 2006; Meneveau 2012; Meyers& Meneveau 2012; Yang, Kang & Sotiropoulos 2012; Abkar & Porté-Agel 2013; Stevens,
Gayme & Meneveau 2016) for more information. Despite the great advantage of thesemodels in capturing the farm–atmosphere interaction, the concept of an infinite wind farm
can be only regarded as an asymptotic case that resembles what very large wind farms may
tend to approach. More importantly, these models fail to offer any insight into the wakeof the wind farm due to their core assumption that the wind farm extends infinitely in thestreamwise direction.
The other group of existing models, which we classify within the broad category ofmulti-scale models, strive to leverage the benefits of both large-scale farm and small-scale
single-turbine modelling. Within this category, different approaches have been adopted
to model wind farm flows. Stevens, Gayme & Meneveau (2015) and subsequent works(e.g. Shapiro et al. 2019; Starke et al. 2021) coupled the infinite-farm approach with thesingle-turbine approach by matching the predicted mean velocity at the turbine hub height.
Other studies coupled the wind farm scale with the turbine scale through a parameter
called the farm induction factor (e.g. Nishino & Dunstan 2020; Kirby, Nishino & Dunstan2022) with more recent works modelling blockage effects as well (e.g. Legris et al.
2022). In another type of multi-scale model, the exchange of energy between the layer
consisting of wind turbines and the overlaying boundary layer was parametrised using theclassical entrainment theory (e.g. Luzzatto-Fegiz & Colm-cille 2018; Bempedelis, Laizet& Deskos 2023). Other multi-scale models characterised the farm–atmosphere interaction
and farm-scale blockage effects caused by meso-scale phenomena such as gravity waves(Allaerts & Meyers 2019; Stipa et al. 2023). The coupling between the different scales inthese models usually involves an iterative process or the numerical solution of governing
equations. Moreover, the focus of the majority of the models discussed above is farm power
production or the flow field within the wind farm, and less attention has been paid to thewake of the entire farm.
In this work, we propose a new category of wake models by considering a semi-infinite
wind farm, i.e. a wind farm that extends infinitely in the lateral direction but has a finite
size in the streamwise direction. The infinite lateral extent of the wind farm allows us toperform lateral averaging, which significantly simplifies the flow’s governing equations
and leads to a closed-form explicit solution without the need for using an iterative
approach. The finite length of the wind farm also makes it possible to systematicallymodel the wake of the entire farm. A schematic of the semi-infinite farm modelling incomparison with single-turbine modelling (i.e. finite approach) as well as infinite-farmmodelling is shown in figure 1. A particular focus of this work is given to the predictionof the deflection of the farm wake. Predicting the magnitude of the farm wake deficit isimportant but not sufficient. The wake deflection also needs to be quantified to determine
whether the wake of a wind farm may impinge on a downstream farm. In general, the985 A43-3

M. Bastankhah and others(b)(a)
FiniteInfinite- Increased aerodynamic roughness- Top-down energy transport- Single-wake modelling- Wake superposition technique

ln(z)
ug

ug
Vertical momentum
transportx(c)Semi-infinite- Infinite width, finite length- Predict streamwise and spanwise deficit

Farm regionyWake region

sy
sx

x

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 1. Different approaches used to model wind farm flows. (a) Modelling each individual wind turbine
wake and then using superposition techniques to account for cumulative wake effects. (b) Modelling a windfarm that is extended to infinity in both streamwise and spanwise directions as an added aerodynamic surface
roughness. (c) Modelling a wind farm that is extended to infinity in the lateral direction but has a finite length
in the streamwise direction.

wake deflection is mainly caused by (i) meso-scale phenomena such as Coriolis force (andits by-product wind veer) and (ii) yaw misalignment. The latter has recently received agreat deal of attention because of its importance in wake steering strategies (e.g. Fleming
et al. 2017; Bastankhah & Porté-Agel 2019; Howland, Lele & Dabiri 2019; Campagnolo
et al. 2020). The former, however, is mostly overlooked in prior modelling works. Whilethe deflection of a single-turbine wake due to Coriolis force is expected to be negligible(Mohammadi et al. 2022), several studies mainly based on numerical simulations haveunderpinned the importance of the farm wake deflection caused by Coriolis force (e.g.
van der Laan et al. 2015; Abkar, Sharifi & Porté-Agel 2016; Allaerts & Meyers 2017;
Eriksson et al. 2019; Gadde & Stevens 2019). Interestingly, there has not been a universal
agreement in the literature with regards to the direction of the wake deflection causedby the Coriolis force. van der Laan & Sørensen (2017) argued that this is due to the factthat the Coriolis force has two effects on the wake deflection. The direct effect turns thewake in the anticlockwise direction (seen from top) in the northern hemisphere, while theindirect effect (through wind veer) rotates the wake in the clockwise direction. Gadde &Stevens (2019) explained this phenomenon based on the direction of the vertical turbulent
fluxes in the entrance region in comparison with those in the wake. One of our objectives
with this new modelling framework is to capture the conflicting influence of the Coriolis
force on the farm wake. This is achieved by concurrently solving momentum equations inboth streamwise and spanwise directions. The outcome is a simple one-dimensional modelthat predicts both laterally averaged streamwise and spanwise velocities at the turbine hubheight within and downwind of the wind farm.
The rest of the paper is organised as follows. Section 2 develops the laterally averaged
Reynolds-averaged Navier–Stokes equations. Section 3 describes the high-fidelity
numerical simulations used in this study to validate the developed model. Section 4
985 A43-4

Fast-running physics-based wake model for a wind farm
discusses the budget analysis that is conducted to identify dominant terms in the
momentum equations. The farm wake model is then developed in § 5. Results are discussed
in § 6, and finally a summary is provided in § 7.
2. Streamwise and spanwise laterally averaged Reynolds-averaged Navier–Stokes
equations

We start by writing the steady-state Reynolds-averaged Navier–Stokes equation for
high-Reynolds-number flows (i.e. negligible friction forces) using Einstein notation. For
simplicity, we non-dimensionalise all variables and equations throughout this paper using
a selection of scales based on the incoming flow and turbine characteristics. All spatial
dimensions are normalised by the turbine rotor diameter D. All velocities ui are normalisedby the incoming velocity Uh at the turbine hub height zh . Static pressure p is normalisedby ρ Uh2 , where ρ is the air density. The dimensionless Coriolis frequency fc is defined as
Dfc= 2Ω sin φ,
Uh

(2.1)

where Ω = 7.2921 × 10−5 rad s−1 is the rotation rate of the Earth and φ is the latitude.
Note that the dimensionless Coriolis frequency fc defined in (2.1) represents the ratioof Coriolis force to inertial force. This dimensionless parameter is in fact the inverse of
the Rossby number Ro that is commonly used in geophysical studies concerning flows inoceans and atmosphere (e.g. van der Laan et al. 2020).
The dimensionless form of the Reynolds-averaged Navier–Stokes equation reads as(Stull 2009) 

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press∂ui uj
∂ p̄
∂ ūi
ūj
= εij3 fc ūj −−− f̄i ,
∂xj
∂xi
∂xj

(2.2)

where ui is the velocity component in the xi direction with i = 1, 2, 3 corresponding to
the streamwise x, spanwise y and vertical z directions, respectively. Overbar denotes time
averaging, and the turbulent fluctuating velocity is ui = ui − ūi . The permutation symbol
is denoted by εijk . Moreover, f̄i is the time-averaged component of the turbine forces perunit volume, non-dimensionalised by ρ Uh2 /D. Note that gravitational forces are neglectedin (2.2) as only the streamwise and spanwise momentum directions are of interest in thisstudy.
We separate the static pressure p in (2.2) into that due to the background driving pressure
gradient pg and that due to the presence of turbines pt . The former is dictated by the force
balance in the geostrophic layer on top of the ABL, while the latter is due to the pressure
drop across the rotor disc and its recovery to the free-stream pressure downstream. The
background driving pressure gradient pg can be written in terms of the geostrophic windspeed ūg (Stull 2009). This simplifies (2.2) to ∂ui uj∂ ūi∂ p̄t= −εij3 fc (ūgj − ūj ) −−− f̄i .
ūj
∂xj
∂xi
∂xj

(2.3)

Different forms of spatial averaging such as surface or volumetric averaging have beenfrequently performed in prior studies on vegetation canopies or infinite wind farms (e.g.
Calaf, Meneveau & Meyers 2010; Moltchanov, Bohbot-Raviv & Shavit 2011; Bai, Katz &Meneveau 2015; Goit & Meyers 2015). In this study, however, we perform lateral averaging(averaging only along the y direction), because our aim is to determine how flow quantities
985 A43-5

M. Bastankhah and othersat the hub-height level evolve along the streamwise direction x. Mathematically, this may
be defined for an arbitrary variable ψ as follows: L
1ψ̄(x, z) = lim
ψ̄(x, y, z) dy,(2.4)
L→∞ 2L −L
where  indicates lateral averaging and [−L, L] is the lateral range over which theaveraging is performed. The lateral fluctuation is defined as ψ  = ψ − ψ and bydefinition ψ   = 0. By performing the lateral average on (2.3), we obtain ∂ui 
∂pt  ∂ui uj  ∂ui  uj  uj 
= − εij3 fc (ugj  − uj ) −−−−  f̄i  .

∂xj
∂xi
∂xj
∂xj



          
T
A

C

P

R

(2.5)

D

The terms outlined in (2.5) are
A Advection of momentum by mean flow.
C Coriolis term.
P Pressure gradient due to the presence of wind turbines.
R Reynolds stress gradients.
D Dispersive stress gradients.
T Turbine forcing.
These terms are discussed in more detail in § 4. The dispersive stress term ui  uj   that isthe product of spatial fluctuations in the lateral direction arises in (2.5) as a result of lateral
averaging. It is also noteworthy that because of lateral averaging any terms including∂/∂xj in (2.5) must be zero if j = 2 (i.e. ∂/∂y = 0).

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

3. Numerical set-up: LES

Large-eddy simulations were performed using the open source software OpenFOAM(version 2.3.1) in conjunction with the Simulator fOr Wind Farm Applications (SOWFA)
project libraries (Churchfield et al. 2012a) developed by the US National Renewable
Energy Laboratory (NREL). The atmospheric solver used in SOWFA is called ABLSolver,which is a transient solver for turbulent flows of incompressible fluids and considers theBoussinesq approximation for buoyancy effects (Churchfield et al. 2012a).
A precursor–successor approach has been utilised to develop a conventionally neutral
ABL flow for the simulation. The Coriolis force is calculated for φ = 55.52◦ , which is a
representative value for a wind farm in the North Sea (Hansen et al. 2012). In SOWFA,
a prescribed streamwise velocity (Uh = 8 m s−1 ) and wind direction (ϕ = 270◦ ) at the
turbine hub-height level can be achieved by adjusting the magnitude and direction of thedriving pressure gradient. A capping inversion with a lapse rate of 0.05 K m−1 is imposedat the top of the boundary layer covering the heights from 700 to 800 m. The height at thebottom of the capping inversion, denoted by H, is defined as the thickness of the ABL.
The geostrophic layer above the capping inversion has a lapse rate of 0.003 K m−1 . The
inclusion of the capping inversion helps to slow the vertical growth of the boundary layer
with time in neutral conditions (Churchfield et al. 2012b). Due to the assumption of the
fixed height of the capping inversion, however, the vertical displacement of the flow above
the farm does not generate gravity waves in the capping inversion. This simplification may
lead to errors in cases where gravity waves induce non-negligible pressure gradients at
985 A43-6

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Fast-running physics-based wake model for a wind farm
the turbine hub-height level (Allaerts & Meyers 2017, 2019; Lanzilao & Meyers 2023;
Stipa et al. 2023). The precursor simulations are run without the turbines for a periodof 10 h (36 000 s) to obtain a quasi-steady state. Next, the inlet conditions are recorded
for a period of 9000 s to be fed into the successor simulation with turbines. The farm
flow statistics are calculated for the last hour of the simulations. The convective terms are
discretised using a second-order central difference scheme for the precursor and a local
blend between linear (second-order) and upwind (first-order) schemes for the successor
simulation depending on the cell size. This scheme uses 80 % linear and 20 % upwind
in proximity of the turbines and 100 % linear in the rest of the domain. For temporal
discretisation, an implicit second-order backward scheme is used. For the diffusion term,a Gauss linear second-order scheme is implemented using a non-orthogonality correction
for surface-normal gradients. Subgrid-scale stresses are modelled using a one-equation
turbulent-viscosity model (Yoshizawa 1986).
The turbine used in this study is NREL-5MW with a hub height of 90 m and a rotor
diameter D of 126 m (Jonkman et al. 2009). The turbines are modelled as an actuator discwith no rotation and a constant thrust coefficient of 0.776. This value of CT was foundbased on blade element momentum simulations of a turbine rotor with Uh of 8 m s−1
and tip-speed ratio of 7.55 (Navarro Diaz et al. 2022). The body forces are spread acrossthe rotor plane uniformly as axial forces. The equivalent inflow velocity is unknown for
turbines that are subject to the wakes of upstream turbines, so a calibration table is used torelate the average velocity on the disc with the unperturbed inflow velocity (van der Laanet al. 2015).
The same domain size is used for both precursor and successor simulations. It extends
1000 m in the vertical direction. The first row of turbines is placed 15 rotor diameters(1890 m) downstream of the inlet, and the domain extends for 357 turbine diameters(approximately 45 km) after the last turbine row. A schematic of the computational domain
is presented in figure 2. It is worth noting that the distance between the inlet and the firstrow of turbines is relatively short in these simulations. Our aim is to maximise the useof our computational resources to capture a very large extent of the farm wake. However,this may lead to an underestimation of the velocity slowdown in the upwind region caused
by farm-scale blockage effects, as discussed in the recent parametric study by Lanzilao &Meyers (2023).
In the precursor simulations, grid cells are 21 m (i.e. D/6) long in the streamwise andlateral directions, and their height in the vertical direction grows with distance from theground (from 2.5 to 60 m at the top of the domain). In the successor simulations including
wind turbines, the mesh is refined in two steps. Each refinement halves the cell size. First,
in a zone containing the wind farm and its downstream region, the mesh size is reduced
to 10.5 m (i.e. D/12) in the streamwise and lateral directions. This refined region starts
1260 m (i.e. 10D) upstream of the first turbine row to ensure that eddy structures are fully
developed in the new refined mesh before reaching the wind turbines. In close proximity
to the wind turbines, the mesh is further refined by a factor of two (i.e. cell size of 5.25 m
or D/24 in the streamwise and lateral directions) to capture strong velocity gradients inthis region.
Precursor simulations use cyclic boundary conditions at the inlet, outlet and sides. The
nearest turbines to the sides are placed such that it resembles the infinite extent of thewind farm in the lateral direction. For instance, with a lateral spacing of 4D, there is a 2D
distance from each side as shown in figure 2. At the ground, a wall boundary condition
with a prescribed roughness length based on the Schumann–Grotzbach formulation is
implemented (Schumann 1975). At the domain top, a slip boundary condition is imposedfor velocity and a fixed gradient for temperature. For successor simulations including wind
985 A43-7

M. Bastankhah and others

357D ≈ 45 km
sy = 4D
sx = 7D

2D

2D

z
y

x

1000 m
15D = 1890 m

Figure 2. Schematic of the computational domain for the A0 case. Turbines are shown by black circles andthe rotor diameter is denoted by D.

Case

Description

Layout

sx

sy

z0,0

N

u∗

I0 (%)ϕ (deg.)

A0
S0
AS
AD
AR

Aligned baseline
Staggered baseline
Aligned short
Aligned dense
Aligned roughAligned
Staggered
AlignedAlignedAligned

7
7
7
5
7

4
4
4
3
4

1.6 × 10−6
1.6 × 10−6
1.6 × 10−6
1.6 × 10−6
1.6 × 10−4

8
8
4
8
8

0.00370.00370.00370.00370.0059

5.6
5.6
5.6
5.6
8.2

−3
−3
−3
−3
−3

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Table 1. Summary of LES for different semi-infinite wind farms. Here sx and sy are, respectively, streamwise
and spanwise inter-turbine spacing, normalised by the rotor diameter D. The surface roughness normalised by
D is shown by z0,0 . The number of turbine rows is denoted by N and u∗ is the friction velocity normalised by
Uh . The incoming turbulence intensity at the hub height is shown by I0 and ϕ is the change in the incoming
wind direction across the turbine rotor (i.e. from the bottom-tip height to the top-tip height).

turbines, the inlet uses the data from the precursor while a zero-gradient condition isapplied at the outlet. The domain’s sides, lower and upper parts have cyclic, wall andslip boundary conditions, respectively. In total, five simulations were performed to study
the effect of wind farm layout, wind farm length, inter-turbine spacing and the incoming
turbulence level on farm wake flows. The details of these simulations are summarised intable 1. In this table, sx and sy are, respectively, streamwise and spanwise inter-turbine
spacing, normalised by the rotor diameter D. The surface roughness normalised by D
is shown by z0,0 . The number of turbine rows is denoted by N, and u∗ is the friction
velocity normalised by Uh . The incoming turbulence intensity at the hub height is shown
by I0 , and ϕ is the change in the incoming wind direction across the turbine rotor (i.e.
from the bottom-tip height to the top-tip height). As seen in table 1, the first four cases(A0, S0, AS, AD) are subject to a smooth boundary layer with low surface roughness,whereas the incoming boundary layer in the AR case has a higher surface roughness.
The two different inflow boundary-layer profiles are shown in figure 3. The instantaneous
streamwise velocity field u for a portion of the Aligned Baseline (A0) case is also shown in
985 A43-8

Fast-running physics-based wake model for a wind farm
(a)

(b)

4.0

(c)

4.0

(d)

4.0

4.0

3.5

3.5

3.5

3.5

3.0

3.0

3.0

3.0

2.5

2.5

2.5

2.5

z 2.0

2.0

2.0

2.0

1.5

1.5

1.5

1.5

1.0

1.0

1.0

1.0

0.5

0.5

0.5

0.5

0

0.6

0.8

1.0

0

1.2

ū

–4 –2

0

2

0

ϕ (deg.)

510I  (%)

150

SmoothRough1

2

uiujh ×10–3

Figure 3. Spanwise-averaged vertical profiles of inflow characteristics obtained from precursor simulations.
(a) The normalised streamwise velocity ū, (b) the wind direction ϕ, (c) the incoming turbulence intensityof streamwise turbulent fluctuations, and (d) the horizontal
I = σu /Uh , where σu is the standard deviation


turbulent shear stress defined as ui uj h = u w 2 + u v  2 . Horizontal dashed and dotted lines respectively
indicate the turbine hub height and vertical positions of top/bottom blade tips.

u
3.9 × 10–2 0.2

0.4

0.6

0.8

1.0

1.2

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Presszyx

Figure 4. Contours of instantaneous normalised streamwise velocity u for the Aligned Baseline (A0) case.
Turbines are shown by black circles.

figure 4, where the highly turbulent nature of the atmospheric flow and low-speed wakes
are clearly visible.
4. Momentum budget analysis

In this section, the LES data for the Aligned Baseline (A0) case are employed toperform budget analysis on the laterally averaged momentum equations (2.5). This analysis
determines dominant terms in the momentum equations both within and downwind of thewind farm. This serves as a basis for the development of the physics-based model later in
985 A43-9

M. Bastankhah and others

Momentum budget terms

0.015
A1x
A2x
Cx
D1x
D2x

0.010
0.005

R1x
R2x
Px
Residual
D1x + D2x + R1x + Px

0–0.005–0.010–0.015–0.020

0

50

100

x

150

Figure 5. Non-dimensionalised streamwise momentum (4.1) budget at hub height for the Aligned Baseline(A0) case. Locations of turbine rows are denoted by vertical dotted lines. All variables are non-dimensionalisedusing a selection of Uh and D.

§ 5. Note that viscous terms in the momentum equations were found to be multiple ordersof magnitude smaller than any other term, so they are neglected in this analysis. Moreover,regions immediately downstream and upstream of the turbine rows are removed due tosteep flow gradients in these regions. This analysis investigates all terms in (2.5), except
for the turbine forcing which is only relevant at the rotor disc.
4.1. Streamwise momentum equation
Writing (2.5) in the streamwise direction (i.e. i = 1) and neglecting turbine forcing gives∂ū∂ū ū ∂ū+ w̄
= − fc (vg  − v̄) −
ū


  ∂x ∂x   ∂z 
 
A1x

Cx

A2x

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press−∂ū w̄ 
 
 ∂z
D2x

D1x

−

∂u u 
 ∂x
 
R1x

−

∂u w  ∂p̄t −,    ∂x
 ∂z
R2x

(4.1)

Px

where {u, v, w} are respectively velocities in the streamwise, spanwise and verticaldirections. The variation of all terms in (4.1) with respect to x is illustrated in figure 5
until 150 rotor diameters downstream of the wind farm, beyond which limited change isobserved in the flow quantities. As shown in figure 5 and as expected, the residual term ismostly negligible throughout the entire domain.
First, we start with the dominant streamwise advection term A1x = ū∂ū/∂x
representing the advection of streamwise momentum by the streamwise velocity.
A positive value for A1x means wake recovery/flow acceleration, and vice versa.
Approaching the farm, A1x becomes negative approximately 8D upstream of the farm.
This is explained by the presence of an induction region preceding the farm that is caused
by farm-scale blockage effects (Bleeg et al. 2018). Behind the first row of turbines, weobserve that the flow acceleration is still suppressed by farm-scale blockage effects, andA1x continues to be negative until about 3D, where the maximum velocity deficit occurs(i.e. A1x = 0). The induction entrance region of the wind farm can be also illustrated bypositive vertical velocity w̄ shown in figure 6(b).
985 A43-10

Fast-running physics-based wake model for a wind farm
(a)

(b)

ū

v̄, w̄

1.0
0.9
0.8
0.7

0

50

100

150

x

200

250

300

v̄

0.010
0.005
0
–0.005–0.010
0

50

100

150

200

w̄

250

300

x

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 6. Variation of laterally averaged (a) streamwise ū and (b) spanwise v̄ and vertical w̄ velocities
with x. All variables are non-dimensionalised using a selection of Uh and D. The locations of turbine rows aredenoted by vertical dotted lines.

After the second turbine row, A1x becomes immediately positive indicating that the
maximum velocity deficit occurs much closer to the turbine, which is followed by flow
acceleration (i.e. wake recovery). By inspection, the profile of vertical Reynolds stressgradient R2x = ∂u w /∂z follows the profile of A1x, confirming that R2x acts to replenish
wake momentum. In other words, peak flow acceleration (i.e. maximum A1x) is observed,
approximately where vertical momentum transport due to turbulence is also maximum.
Note that terms on the right-hand side of (4.1) that are negative in figure 6 promote wakerecovery and vice versa. The greater proportions of turbulent vertical momentum transport
in later rows is also evident in figure 5, occurring due to increased flow shear from greater
velocity deficits, as observed in figure 6(a). It is worth reminding ourselves that according
to (4.1), the gradient of the Reynolds stress is responsible for wake recovery, as opposed tothe common assumption that the mere presence of Reynolds stress promotes wake recovery(van der Laan, Baungaard & Kelly 2022).
Figure 5 also shows that the other advection term A2x = w̄∂ū/∂z has minimal impact
on momentum transport within the domain. Moreover, although not discernible in figure 5,
the Coriolis term Cx = fc (vg  − v̄) is negative across the domain as expected in thenorthern hemisphere, but like A2x, its value is negligible compared with the dominant
terms.
The normal Reynolds stress gradient R1x = ∂u u /∂x is illustrative of the rate ofchange of turbulence level (intensity) with x. Term R1x in figure 5 indicates the turbulence
level increases behind each turbine row, peaking about 3D–5D downstream (i.e. whereR1x = 0), before decreasing on the approach to the next row. This is in agreementwith prior studies observing peak turbulence intensity occurring a few rotor diameters
downstream of an individual turbine (Wu & Porté-Agel 2012). In the wake of the farm,
R1x quickly approaches zero as the turbulence decays to its background level.
Figure 5 shows that the turbine pressure gradient Px = ∂p̄t /∂x is significant within
the entire farm. From actuator disc theory (Manwell, McGowan & Rogers 2010), we knowthat after the pressure increase upwind of turbines, there is a sudden pressure drop as
the turbine extracts energy from the flow (not shown in figure 5). This is followed by apressure increase as wake recovery occurs. Figure 5 shows the fast recovery of pressure
downwind of each turbine row indicated by positive P1x. The value of P1x (i.e. rate ofpressure increase) decays with x until it increases again due to the induction region ofsubsequent rows. It is worth noting that the variation of pressure is often neglected in wakemodels, but this figure and other recent studies (e.g. Bastankhah et al. 2021) highlights
the importance of this term. According to figure 5, within the wind farm, this term iseven comparable to other dominant terms (e.g. A1x and R2x) in the momentum equation.
985 A43-11

M. Bastankhah and others

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

The term Px, however, decays quickly in the wake of the wind farm as the pressure
approaches its free-stream value.
Some of the most significant terms in figure 5 are the dispersive stress terms which
are the product of deviations from the lateral averaging, and described as the tortuous
streamlines induced by flow obstacles (Moltchanov et al. 2011). Dispersive stresses
are correlated to obstacle density. Sparsely populated obstacle fields display increased
dispersive stresses, due to greater disparity between flow over the obstacles and the mean
flow (Moltchanov et al. 2011). For instance, dispersive stresses are expected to be greater inaligned wind farms (shown in figure 5) than in staggered ones (not shown here) although
D1x may be still considerable within a staggered wind farm. The term D1x represents
the amount of streamwise momentum transport caused by flow inhomogeneity. Moreprecisely, it represents the rate of change in inhomogeneity with respect to the laterally
averaged flow. According to figure 5, D1x is negative within the wind farm indicating thehomogeneity of the flow is increasing. This occurs as the wake recovers due to vertical
momentum transport into the wake, reducing the magnitude of the spatial fluctuations ofstreamwise velocity (i.e. ū ). Accordingly, the location of minimum D1x (i.e. maximum
rate of approaching homogeneity) is correlated with where maximum vertical momentum
transport R2x and maximum wake recovery rate A1x occur. As mixing increases so does
the homogeneity, causing the reduction of the magnitude of D1x displayed in figure 5.
Despite its importance within the farm, D1x decays sharply in the wake of the wind farm,where individual turbine wakes merge and form a holistic farm wake. This is discussed in
more detail in § 6. Finally, we investigate the variation of D2x = ∂ū w̄ /∂z. This term
essentially quantifies the vertical transfer of streamwise momentum caused directly by the
non-uniformity of the time-averaged wind farm flow field. Figure 5 shows that this termis clearly smaller than the other dispersive term D1x. It is also interesting to note that thisterm is mainly positive within the wind farm. This indicates that D2x acts against wakerecovery, in contrast to its turbulent counterpart R2x.
4.2. Spanwise momentum equation
Even though the terms of the spanwise momentum equation are of smaller magnitude
than their streamwise counterparts, they are examined here due to their importance to thewake’s trajectory. Writing (2.5) in the spanwise direction (i.e. i = 2) yields∂v̄∂v̄  ū  ∂v̄  w̄ ∂v̄+ w̄
= − fc (ū − ug ) −
−
ū

  ∂x ∂x   ∂z    ∂z
 
A1y

Cy

A2y

−

D1y

∂v  u  ∂v  w −.
 ∂x   ∂z
 
R1y

D2y

(4.2)

R2y

Variations of terms in (4.2) with the streamwise distance x are shown in figure 7. Dueto their small values, the terms in (4.2) are more prone to numerical errors, which may
explain why their variations are more oscillatory and less smooth compared with theirstreamwise counterparts.
First, we start with the Coriolis term Cy = fc (ū − ug ). The dominance of thisterm varies across the domain. From figure 7, Cy is dominantly negative within thefarm due to the streamwise flow deceleration, which according to (4.2) causes positive
A1y = ū∂v̄/∂x and thereby positive v̄, as observed in figure 6(b). This deflects the985 A43-12

Fast-running physics-based wake model for a wind farm(×10–3)
A1y
A2y
Cy
D1y

Momentum budget terms

1.0

0.5

D2y
R1y
R2y
Residual

0

–0.5

–1.0
0

50

100

150

200

250

300

x
Figure 7. Non-dimensionalised spanwise momentum (4.2) budget at hub height for the Aligned Baseline (A0)case. Locations of turbine rows are denoted by vertical dotted lines. All variables are non-dimensionalisedusing a selection of Uh and D.

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

wake to the left, which can be described as an anticlockwise flow rotation viewed fromthe top. In the wind farm wake, figure 7, however, shows that with flow acceleration
the strength of the Coriolis force decays. This is where R2y = ∂v  w /∂z that is thevertical turbulent entrainment of veered momentum from above becomes more dominant.
Consequently, this changes the sign of the advection term A1y to negative as shown infigure 7, and therefore the far wake starts deflecting to the right (i.e. clockwise wakerotation viewed from top). In other words, the term A2y turns the wind at the hub height
towards the wind direction at higher altitudes. This ultimately leads to a negative spanwise
velocity as depicted in figure 6(b). This interesting phenomenon is discussed in more detailin § 6. It is also noteworthy that the magnitudes of the second advection term A2y, thedispersive stress terms D1y and D2y and also R1y are small, especially in the farm wakeand can be neglected.

4.3. Approximate form of momentum equations
From the analysis in § 4.1, amongst others the dispersive stress (D1x) and pressure(Px) terms are evidently not negligible in the streamwise momentum equation, at least
within the farm region. However, despite the evident importance of these two terms, thesummation of the four terms D1x, D2x, R1x and Px – which are challenging to model – israther small. This is illustrated by the dashed light green colour in figure 5. The combined
value of these terms is negligible in the wake of the farm. Within the farm, the combined
value is not negligible but smaller than the individual dispersive D1x and pressure Px
terms. The term D1x is negative, while Px is positive; therefore, to some extent, they
cancel each other out. For simplicity, we thus omit these terms from our model developedin § 5. Moreover, as discussed in § 4.2, it is apparent that in the spanwise direction thedominant terms are (i) the spanwise momentum advection by the streamwise velocity A1y,(ii) the Coriolis term Cy and (iii) the vertical turbulent transport of veering wind R2y.
Therefore, the approximate forms of the momentum equations including turbine forcing
985 A43-13

M. Bastankhah and others
can be written as
ū

∂ū
∂u w ≈ −fc (vg  − v̄) −−  f̄x ,∂x∂z(4.3a)ū∂v̄
∂v  w 
≈ +fc (ug  − ū) −−  f̄y .
∂x∂z(4.3b)In the following section, we simplify (4.3) to develop a system of ordinary differential
equations that can be solved mathematically for a semi-infinite wind farm.
5. Derivation of physics-based fast-running farm wake model

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

5.1. Definition of semi-infinite wind farm
We start by assuming a semi-infinite wind farm which is infinite in the lateral direction but
has a finite length in the streamwise direction. A right-handed Cartesian coordinate system{x, y, z} aligned with the incoming wind at the turbine hub height zh is adopted such thatx is in the direction of the incoming wind, y represents the horizontal direction normal tox and z measures the height from the ground. The total number of wind turbine rows isdenoted by N. Wind turbines in the nth row are abbreviated to WTn s, where the subscript
n = {1, 2, . . . , N} shows the row number labelled based on the streamwise position (i.e.
ranging from n = 1 for the first row to n = N for the last row). The WTn s are assumed tohave the same values of thrust coefficient CT,n and yaw angle γn , which may, however, bedifferent from CT,m and γm if n =/ m. In other words, turbines in different rows may havedifferent operating conditions. While the lateral spacing sy between turbines is assumed
to be the same for all rows, the streamwise spacing between consecutive rows may bevariable. The arbitrary streamwise positions of turbine rows is quite advantageous, as for
instance the developed model can be even applied to a cluster of wind farms at once (notdone in this study). Furthermore, turbine rows may have lateral offset with respect to eachother as shown in figure 8. The lateral position of WTn s is yn + ksy , where k = −∞ . . . ∞.
5.2. Turbine force
The normalised aerodynamic force f̄ exerted by wind turbines is given by
N 
∞

1
f̄x = −
CT,n ū2h,n cos(γn )δ(x − xn )H(0.52 − [( y − yn − ksy )2 + (z − zh )2 ]),2
n=1 k=−∞(5.1a)f̄y = −N ∞
1CT,n ū2h,n sin(γn )δ(x − xn )H(0.52 − [( y − yn − ksy )2 + (z − zh )2 ]),2
n=1 k=−∞(5.1b)
where ūh,n is the local incoming hub-height velocity for WTn s. In other words, it is thetime-averaged velocity at the location of WTn s’ rotor centre (i.e. (x, y, z) = (xn , yn , zh )) in
the absence of WTn s. Here δ(x) is the Dirac delta function and H(x) is the Heaviside step
function, which is defined as H(x) = 0 for x ≤ 0 and H(x) = 1 for x > 0. We then apply
the lateral averaging discussed in § 2 to obtain laterally averaged turbine forces at the hub
985 A43-14

Fast-running physics-based wake model for a wind farm
WT3
y

WT2

WT1

WT4

sy
sy

sy
sy
y4
y3
y2
y1

x1

x2

x3

x4

x

Figure 8. Schematic of a semi-infinite wind farm. Turbines in row n are denoted by WTn . The lateral spacing
sy between turbines is assumed to be constant for the whole wind farm. However, the streamwise spacing canbe variable, and moreover turbine rows may be laterally shifted with respect to each other.

height z = zh as follows:
1 
CT,n ū2h,n cos(γn )δ(x − xn ),2sy(5.2a)1  f̄y  = −CT,n ū2h,n sin(γn )δ(x − xn ).
2sy(5.2b)N f̄x  = −n=1
N

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

n=1

5.3. Simplifying Reynolds shear stress terms in (4.3)
The Boussinesq eddy-viscosity hypothesis has been used in previous studies (Belcher,
Jerram & Hunt 2003; Luzzatto-Fegiz & Colm-cille 2018) to model spatially averaged
Reynolds stresses based on spatially averaged velocity gradients:
⎡
⎤
⎢
⎢

1
ui uj  = −2νt ⎢⎢⎣2⎥∂ui  ∂uj  ⎥
⎥
+
∂xj
∂xi ⎥

⎦

(for i =
/ j),

(5.3)

Sij

where Sij is the dimensionless rate of strain tensor and νt is the turbulent viscosity
non-dimensionalised by Uh D. For i = 1 and j =/ 1, ∂ui /∂xj is an order of magnitude larger
than ∂uj /∂xi (Tennekes & Lumley 1972), so one can write∂u w ∂v  w 
∂ 2 ū∂ 2 v̄
≈ −νt
≈
−ν
,
.
(5.4a,b)
t
∂z
∂z
∂z2
∂z2
We first assess the validity of the turbulent-viscosity hypothesis using our LES data. To doso, at a given streamwise position, we examine whether −u w  is linearly proportional to985 A43-15

M. Bastankhah and others
x = 2.5sx–uw

(a)x = 6.5sx(b)x = 1.5xN(c)x = 2xN

(d)

0.004

0.004

0.004

0.004

0.002

0.002

0.002

0.002

0

0

0

0

4
3

0

0.1

∂ū/∂z

0.2

0

0.1

0.2

0

∂ū/∂z0.1∂ū/∂z

0.2

2

z

1
0

0.1

0.2

∂ū/∂zhttps://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 9. Distribution of −u w  against ∂ū/∂z at different heights for four different streamwise positionsin the A0 case, where (a,b) are in the farm region and (c,d) are in the wake region, and xN = 7sx is the farm
length.

∂ū/∂z, according to (5.3). Figure 9 shows values of −u w  and ∂ū/∂z at different
heights and streamwise positions. Results are shown for the A0 case, but they look
qualitatively similar in other cases (not shown here). Data are plotted for a range ofz = 0.25 (shown in light blue) to z = 4.57 (shown in magenta). Results in figure 9 areshown for four streamwise locations, with the first two in the farm region and the othertwo in the wake region: between WT3 and WT4 (x = 2.5sx ) (figure 9a), between WT7 and
WT8 (x = 6.5sx ) (figure 9b), half of the farm length downstream in the wake (x = 1.5xN ),where the farm length is xN (figure 9c) and an entire farm length downstream in thewake (x = 2xN ) (figure 9d). Lines are fitted to the data at heights above the hub height.
Figure 9(a,b) shows that within the farm region, the eddy-viscosity assumption seems to
be a valid approximation for the entire vertical domain plotted in the figure, except for theregion close to the ground. In the farm wake, however, this assumption is only valid at
upper heights as shown in figure 9(c,d). Given that the height at which the eddy-viscosity
assumption appears reasonable in the farm wake seems to grow with downstream distance,
one possible explanation for this could involve the development of the secondary internalboundary layer (IBL) downwind of the wind farm due to the transition from rough tosmooth terrain (see § 5.5 for more discussion on IBLs). However, further research is
required to fully elucidate this phenomenon. Nevertheless, even in the farm wake, theturbulent-viscosity hypothesis is still valid at upper heights, where the top-down transport
of energy by Reynolds shear stresses occurs, so we use the turbulent-viscosity assumption
in this work to simplify the laterally averaged momentum equations.
The turbulent viscosity νt is then decomposed into two parts: one that is due to theambient atmospheric turbulence denoted by νt,0 and one shown by νt,f corresponding to
the turbulence added by the wind farm, and it varies by x. In other words,
νt (x) = νt,0 + νt,f (x).

(5.5)

Specifying values of the ambient turbulent viscosity νt,0 and the farm turbulent viscosity
νt,f is deferred to § 5.5.
5.4. Mathematical solution of velocity deficit equations
Now, we develop and solve equations for the variation of laterally averaged velocity deficit
in both streamwise Ud and spanwise Vd directions, defined as
Ud (x) = U0 − U(x),(5.6)

Vd (x) = V0 − V(x).

(5.7)

In the above equation and hereafter, U(x) = ū(x, z = zh ), V(x) = v̄(x, z = zh ). The
subscript 0 denotes the flow in the absence of the whole wind farm. Let us recall that
985 A43-16

Fast-running physics-based wake model for a wind farm
the coordinate system is defined based on the incoming wind direction at the hub height(see § 5.1) and all velocities in this work are non-dimensionalised by Uh , so U0 = 1 and
V0 = 0. The latter means that Vd = −V in this work. For generality, however, we write
equations for Vd so the model can still be used for a different coordinate system where
V0 =
/ 0.
The magnitude of the local velocity deficit caused by each individual turbine can besignificant especially in the near-wake region (e.g. Zhang, Markfort & Porté-Agel 2012;
Bastankhah & Porté-Agel 2017). However, laterally averaged values of the velocity deficit
are fairly small, as is later shown in § 6. Therefore, we linearise (4.3) by replacingū in the advection term with the incoming velocity U0 = 1. Moreover, using the
turbulent-viscosity hypothesis discussed in § 5.3 to simplify (4.3), we obtain∂ 2U∂U≈ −fc (Vg − V) + νt 2 −  f̄x ,∂x∂z(5.8a)∂ 2V∂V≈ +fc (Ug − U) + νt 2 −  f̄y .
∂x∂z(5.8b)In the absence of the wind farm, (5.8) is simplified to0 ≈ −fc (Vg − V0 ) + νt,0

∂ 2 U0,∂z2(5.9a)∂ 2 V0.
∂z2(5.9b)

0 ≈ +fc (Ug − U0 ) + νt,0

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

It is worth noting that (5.9) was solved in Ekman (1905) for different altitudes to
describe the well-known Ekman spiral. If we subtract (5.8) from (5.9) and also use thedimensional analysis of ∂ 2 Ud /∂z2 ∝ Ud /l2 , where l ≈ 1 (i.e. length scale comparable torotor diameter), we obtaindUd≈ + fc Vd − c1 νt Ud +dx   
Coriolis

recoverydVd≈ − fc Ud − c1 νt Vd +dx   
Coriolis

recovery f̄x 

+ Cx ,
turbine forcing

shear

 f̄y 

+ Cy ,
turbine forcing(5.10a)(5.10b)veer

where c1 is a constant. The term Cx = −νt,f ∂ 2 U0 /∂z2 is related to the rate of incoming
shear, so it is called shear in (5.10). Likewise, Cy = −νt,f ∂ 2 V0 /∂z2 is related to the rate ofincoming veer, and it is called veer in (5.10). Both Cx and Cy also depend on the amount of
turbulence generated by the wind farm through νt,f . Values of Cx and Cy are determined asa function of atmospheric and farm conditions later in § 5.7, and for now they are assumed
to be known.
If we insert (5.2) into (5.10), the mathematical solution presented below for this systemof ordinary differential equations can be obtained. The solution written in (5.11) is exact
for a constant νt,f . For a well-defined function of νt,f (x) with an arbitrary distribution,
this provides an approximate solution with insignificant error with respect to the exact
985 A43-17

M. Bastankhah and others
numerical solution (not shown here) of the ordinary differential equation system:xCx(1 − e−c1 0 νt dx )H(x) +Ud,n (x),c1 νt(5.11a)xCy(1 − e−c1 0 νt dx )H(x) +Vd,n (x),c1 νt(5.11b)N

Ud (x) ≈n=1
N

Vd (x) ≈n=1

where Ud,n (x) and Vd,n (x) are, respectively, values of the laterally averaged streamwise
and spanwise velocity deficit caused by WTn s at x, and they are given by
Ud,n (x) =xCT,n(1 − ηn Ud (xn ))2 e−c1 xn νt dx cos(γn − fc (x − xn ))H(x − xn ),2sy(5.12a)Vd,n (x) =xCT,n(1 − ηn Ud (xn ))2 e−c1 xn νt dx sin(γn − fc (x − xn ))H(x − xn ).
2sy(5.12b)In (5.12), ūh,n is replaced withhttps://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Pressūh,n = 1 − ηn Ud (xn ),(5.13)

where Ud (xn ) = Ud (x = xn ), and we call η the local deficit coefficient, which is definedas the ratio of the local velocity deficit experienced by wind turbines to the laterally
averaged velocity deficit. This coefficient depends on the farm layout configuration, and
it is determined by (5.26) developed later in § 5.6. To compute the ambient turbulent
viscosity νt,0 and farm turbulent viscosity νt,f , (5.14) and (5.19) developed in § 5.5 arerespectively used. Finally, values of Cx and Cy are estimated based on (5.30) in § 5.7.
Once values of Cx , Cy , νt,0 and νt,f are all known, a forward marching scheme in thestreamwise direction is implemented using (5.11) to find the evolution of Ud and Vd withx. The forward marching scheme is stable and not sensitive to the streamwise resolution,
which was tested (not shown here) for a range of values from 0.01D to 1D. The solution in(5.11) is versatile as it can also be used for different distributions of νt , Cx and Cy that canpotentially be developed in future studies.
5.5. Estimation of turbulent viscosity νt = νt,0 + νt,f
In general, turbulent viscosity can be written as a product of a turbulence velocity scale ûand a turbulence length scale (i.e. mixing length) l̂ (Tennekes & Lumley 1972). Therefore,
to estimate the ambient turbulent viscosity νt,0 and the farm turbulent viscosity νt,f , one
needs to specify suitable values of û and l̂ for each term.
For the ambient turbulent viscosity νt,0 , according to Prandtl’s mixing-length hypothesis(Tennekes & Lumley 1972), we assume û0 = u∗ and l̂0 = κzh , where κ ≈ 0.41 is thevon Kármań constant. Therefore, the ambient turbulent viscosity is given by
νt,0 = κu∗ zh .

(5.14)

It is worth noting that the assumption of l̂0 = κzh is expected to be valid only in the
log-law region of the ABL (Pope 2000). For our simulations the top of the rotor disc is
at a location z/H = 0.22 which is seemingly right at the outer limit of this assumption.
More sophisticated models for the mixing length in ABLs can be used instead (e.g.
985 A43-18

Fast-running physics-based wake model for a wind farm
van der Laan et al. 2020) but for simplicity we retain the simple log-law relationship forl̂0 .
For the farm turbulent viscosity νt,f , we use the height of the IBL that grows above thewind farm as the turbulence length scale. Several studies have shown that the IBL grows
above wind farms following the classical Elliott’s x0.8 power law (e.g. Allaerts & Meyers
2017; Wu & Porté-Agel 2017). According to Wood (1982), the thickness of the IBL, δ, due
to a smooth to rough transition is given byδz0,f= 0.28

0.8

x

,

z0,f

(5.15)

where z0,f is the equivalent roughness length of the wind farm. It is worth noting that the
IBL may become separated from the ground surface as a new IBL starts developing due
to the rough to smooth transition downwind of the wind farm (Oke 1976). For simplicity,
we use the maximum height of the first IBL as the turbulence length scale and do notaccount for the second IBL development. In order to use (5.15), one needs to estimate
the value of z0,f for the wind farm. For an infinite wind farm, several models have been
already proposed in the literature to estimate z0,f (e.g. Frandsen 1992; Calaf et al. 2010;
Yang et al. 2012; Abkar & Porté-Agel 2013). The one suggested by Frandsen (1992) for an
infinite wind farm with uniformly distributed wind turbines states⎛⎞⎜⎜
z0,f = zh exp ⎜
⎜− 
⎝
1

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

2


cft +

κ
κ
ln(zh /z0,0 )⎟⎟⎟2 ⎟ ,⎠(5.16)

where cft = πCT /(4sx sy ), sx is the normalised streamwise spacing between turbine rows
and as a reminder z0,0 is the normalised surface roughness length in the absence of thewind farm. To maintain simplicity, we use the same relationship to estimate z0,f for thesemi-infinite wind farm. To compute cft , we use the average streamwise inter-turbine
N−1spacing for sx (i.e. sx = n=1(x− xn )/(N − 1)), and the average value of thrustN n+1coefficient for CT (i.e. CT = n=1 CT,n /N).
One can use (5.15) in conjunction with (5.16) to estimate the thickness of the IBL over
the wind farm. This relationship is, however, only valid as long as δ is smaller than the
ABL thickness H (Wood 1982). It is known that the IBL growth is capped by the thermal
inversion at the top of the ABL (Oke 1976), especially if the inversion layer has strong
free-atmosphere stratification. In these cases, the inversion layer may act as a ‘lid’ on the
top of the ABL and hinders the growth of the IBL. We therefore artificially limit the growthof the turbulence length scale l̂f using the following relationship:
l̂f =

δ
.
1 + δ/H(5.17)For small values of x, l̂f ≈ δ, while l̂f → H as the value of x → ∞. Variation of l̂f with x
based on (5.17) is shown in figure 10. For simplicity, we here assume that the value of H isconstant. It is, however, important to note that depending on the size of the wind farm andthe level of thermal stratification in the inversion layer, the IBL may lead to the growthof the entire ABL, so in reality H might change with x (Allaerts & Meyers 2017; Wu &Porté-Agel 2017).
985 A43-19

M. Bastankhah and others
1.0

ûf
lˆf
νt, f ∝ ûf lˆf
νt, f (LES – case S0)

0.5
0

0

50

100

x

150

200

250

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 10. Variation of velocity scale ûf , length scale l̂f and farm turbulent viscosity νt,f based on themodelling approached elaborated in § 5.5. Variation of νt,f for the LES data (case S0) is also shown. In the√
figure, ûf is normalised by cft , νt,f is normalised by νt,f (x = x1 + Lf ) and l̂f is normalised by H. Verticaldotted lines show the location of turbine rows, and the vertical dashed line shows x1 + Lf where the turbulent
viscosity νt,f is maximum.

Next, we determine the turbulence velocity scale ûf . Within the wind farm, turbulence
is mainly generated due to the shear caused by turbine forcing. Therefore, inspired by√(Calaf et al. 2010), we use cft U0 , where U0 = 1, to estimate the turbulence velocity
scale ûf within the wind farm. The turbulence generation is expected to peak at somedistances downstream of the wind farm, which is then followed by a decay in turbulence
generation due to wake recovery and reduction of flow shear (Stieren & Stevens 2022).
The generated turbulence in the turbine wake usually peaks at around 5 rotor diameters
downstream (Chamorro & Porté-Agel 2009). The constant turbulence velocity scale ûfcaused by turbine forcing is therefore assumed to be extended from x1 to xN + 5 as shownin figure 10. Further downstream, ûf starts to decay due to the wake recovery. The wakeof a semi-infinite wind farm can be modelled as a two-dimensional wake of a canopy ofroughness elements in a turbulent boundary layer. According to Belcher et al. (2003), thevelocity scale defined based on the maximum velocity deficit in this type of wake flows
decays with x−1 . So in summary, we model the turbulence velocity scale ûf as⎧⎪
0
if x < x1 ,⎪⎪⎨√c
if x1 < x < x1 + Lf ,ft
(5.18)
ûf =
−1
⎪
x
−
x
√
1
⎪
⎪
c
if
x
>
x
+
L
,
⎩ ft
1
f
Lf
where Lf = (xN − x1 ) + 5. Variations of ûf and νt,f are shown in figure 10, where the farmturbulent viscosity νt,f is given by
νt,f = c2 ûf l̂f ,(5.19)

with c2 a constant. As can be seen in figure 10, the farm turbulent viscosity increases
within the farm until it reaches its maximum value five rotor diameters downstream of
the wind farm. It then decreases further downstream until it eventually becomes zero very
far downstream. In general, this is in fairly good agreement with our LES data. As an
example, the variation of νt,f for the Staggered Baseline (S0) case is shown in figure 10.
Fairly similar variations of turbulent viscosity have been also reported in the literature forthe wake of a single wind turbine (Scott et al. 2023).
5.6. Estimation of local deficit coefficient ηnAs discussed earlier, ūh,n = 1 − ηn Ud (x = xn ), where ηn is the local deficit coefficient for
WTn s. A consequence of linearising the momentum equation (5.8) is that wake effects are
985 A43-20

Fast-running physics-based wake model for a wind farm
linearly superposed (Bastankhah et al. 2021). Therefore, one can writen−1ūd,m (xn , yn , zh )
m=1
ηn =
,
n−1


(5.20)

Ud,m (xn )m=1

where ūd,m (x, y, z) is the local velocity deficit at (x, y, z) caused by WTm s and Ud,m (x)is already defined in (5.12). To find ūd,mn , one can use a Gaussian distribution
C exp(−r2 /2σ 2 ) to express the velocity deficit caused by each turbine, where C is
the wake-centre velocity deficit, r is the distance from the wake centre and σ is thecharacteristic wake width. Therefore, for a semi-infinite wind farm, we obtainūd,m (xn ) = lim Cmn
K→∞

K


exp −

k=−K

( yn − ym − ks)2
2
2σmn

,

(5.21)

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

where Cmn is the maximum velocity deficit caused by each WTm at x = xn , k denotes
the column number which varies from −∞ to ∞ and σmn is the width of WTm wakes at
x = xn . Solving (5.21) gives!
"#2√
π( yn − ym )
σmn
σ
mn
ϑ3
, exp −2π2 2,(5.22)ūd,m (xn ) = 2πCmn
sy
sy
sy

l2
where ϑ3 [z, q] is the Jacobi theta function defined as 1 + 2 ∞l=1 q cos(2lz) (Whittaker& Watson 2020). High-level programming languages (e.g. Python) may have a built-infunction to compute the Jacobi theta function (mpmath development team 2023).
The series describing the Jacobi theta function converges rather quickly, so as an
approximation, one can alternatively compute the summation of the first few terms in theseries. It is also noteworthy that ϑ3 [z ± π, q] = ϑ3 [z, q], so yn and ym in (5.22) can be the
spanwise location of any turbines in WTn and WTm , respectively.
By definition Ud,m (x) = ūd,m (x), so we can write
Cmn Ud,m (xn ) = lim
K→∞ Msy
K

k=−K

 +∞

!

( y − ym − ksy )2
exp −
2
2σmn
−∞

"
dy =

√

2πCmn

σmn
.
sy
(5.23)

Therefore, from (5.22) and (5.23), we conclude thatūd,m (xn ) = Ud,m (xn )ϑ3"#!
2π( yn − ym )σmn.
, exp −2π2 2sysy(5.24)

The wake width σmn can be simplified by k∗ (xn − xm ) + σ0 , where the wake expansionk∗ is typically around 0.02–0.04 for offshore conditions and σ0 is around 0.2–0.3(Bastankhah & Porté-Agel 2014). For simplicity, we assume k∗ = 0.025 and σ0 = 0.25 to
985 A43-21

M. Bastankhah and others
reduce (5.24) toūd,m (xn ) ≈ Ud,m (xn )ϑ3"#!
π( yn − ym )(xn − xm + 10)2.
, exp −sy
80s2y

(5.25)

Therefore from (5.20) and (5.25), we conclude
n−1


ηn =

m=1

Ud,m (xn )ϑ3"#!
(xn − xm + 10)2π( yn − ym ), exp −
sy
80s2y
n−1


,

(5.26)

Ud,m (xn )https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

m=1

where n > 1 (for n = 1, η is zero) and Ud,m (x) is defined in (5.12). Figure 11 shows thevariation of the local deficit coefficient ηn with the row number n for both the Aligned
Baseline (A0) and the Staggered Baseline (B0) wind farms. It is interesting to note thatthe value of ηn is always greater than one for an aligned wind farm. This is expected asin this case turbines operate in full-waked conditions. Therefore, the local velocity deficit
experienced by turbines is expected to always be larger than the laterally averaged velocity
deficit. The maximum value of η occurs at the second row where there is the maximum
level of heterogeneity in the farm. Further downstream, due to the wake expansion andflow mixing, the value of ηn decays and approaches a constant value of around two forthis particular wind farm layout. For the staggered wind farm, however, we observe a very
different behaviour. For the second row, the value of ηn is almost zero. This is due to thefact that WT2 s are not affected by the wake of WT1 s due to the staggered layout of the wind
farm. Wake interactions only occur from the third row where there is a sudden increase in
the value of ηn . The value of ηn for the staggered wind farm approaches one for WT4 anddownwind rows. This suggests a fairly uniform distribution of streamwise velocity deep
inside a staggered wind farm. In figure 11(b), we can see how different distributions of thelocal velocity coefficient ηn lead to different distribution of local hub-height velocity ūh,n .
The local hub-height velocity is clearly higher for the staggered wind farm layout whichleads to more power generation. The trend shown in figure 11(b) is similar to the datareported in other works (e.g. Chamorro, Arndt & Sotiropoulos 2011; Stieren & Stevens2022). We have only discussed aligned and staggered layouts here, but an important featureof our developed model is that it is generalisable, through variable local deficit coefficientηn developed in (5.26), to any conceivable wind farm layout provided the layout fulfils therequirements defined in § 5.1.
5.7. Estimation of shear term Cx and veer term Cy
As discussed earlier in (5.10), Cx = −νt,f ∂ 2 U0 /∂z2 and Cy = −νt,f ∂ 2 V0 /∂z2 . Based on(5.9), Cx and Cy can be written asνt,ffc (Vg − V0 ),
νt,0
νt,f
fc (Ug − U0 ).
Cy =
νt,0

Cx = −(5.27a)(5.27b)As a first-order approximation, the vertical changes in the streamwise and spanwise
velocities across the ABL are proportional to G cos θ0 and G sin θ0 , respectively.
985 A43-22

Fast-running physics-based wake model for a wind farm
(a)

4
3

ηn 2
1
0
1

2

3

4

5

6

7

8

(b)
ūh,n = 1 – ηnUd (xn)1.0
Staggered – case S0
0.9

Aligned – case A0

0.8
0.7
0.6

1

2

3

4

5

6

7

8

n
Figure 11. Variation of (a) local deficit coefficient ηn and (b) local hub-height velocity ūh,n with row number
n for the Staggered Baseline (S0) and the Aligned Baseline (A0) cases based on (5.26).

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press
Here, G = Ug2 + Vg2 is the geostrophic wind speed, and the cross-isobar angle θ0 is theangle between the wind direction on the ground surface and the geostrophic wind direction.
One may thus write
Ug − U0 = c3 G cos θ0 ,(5.28a)Vg − V0 = c3 G sin θ0 ,(5.28b)where c3 is a constant. Values of G and θ0 can be obtained from the widely usedgeostrophic drag law, which relates surface properties (e.g. z0 and u∗ ) to geostrophic windspeed on top of the ABL (e.g. Blackadar & Tennekes 1968; Hess & Garratt 2002; van derLaan et al. 2020). For a neutrally stratified ABL, the geostrophic drag law reads as (Liu,
Gadde & Stevens 2021)
A = lnκGu∗−cos θ0 ,z0,0 |fc |
u∗
κG
B=∓
sin θ0 ,u∗(5.29a)(5.29b)where A and B are universal empirical constants, and values of A = 1.8 and B = 4.5 usedin the Wind Atlas Analysis and Application Program (Floors, Troen & Kelly 2018) are
adopted here. The minus sign on the right-hand side of the second equality in (5.29)
relates to the northern hemisphere where fc is positive, whereas the positive sign is forthe southern hemisphere where fc is negative (Liu et al. 2021). From (5.27), (5.28) and
985 A43-23

M. Bastankhah and others
c1

c2

c3

0.85

0.05

0.06

Table 2. Empirical model coefficients used in this study.

(5.29), we can therefore approximate the values of Cx and Cy as follows:
Cx = c3 |fc |Cy = c3 fc

νt,f u∗B,
νt,0 κ


νt,f u∗u∗ln
νt,0 κ
z0,0 |fc |(5.30a)−A .

(5.30b)

While for simplicity the conventional geostrophic drag law relationship is used here
to estimate Cx and Cy , recent works (e.g. Liu & Stevens 2022; Narasimhan, Gayme& Meneveau 2023) that describe the structure of Ekman boundary layer flows can be
implemented in future works.

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

6. Results and discussion

In this section, we compare predictions of the model developed in § 5 with the LES dataelaborated in § 3. The values of model coefficients, namely c1 , c2 and c3 , used in this study
are listed in table 2. It is worth recalling that c1 is the coefficient that is multiplied to νt =
νt,0 + νt,f in the final solution, while c2 is the coefficient within νt,f . An increase of either
c1 or c2 increases the wake recovery rate, but c2 is the one that quantifies the importance of
νt,f over νt,0 . On the other hand, c3 is the coefficient of Cx (i.e. shear term in the streamwise
momentum equation) and Cy (i.e. veer term in the spanwise momentum equation). Given
that Cx is fairly small compared with other terms in the streamwise momentum equation,
the main impact of c3 is on predictions of the spanwise velocity deficit. The reported
coefficients were found manually based on comparing model outputs with the LES data,
but a more systematic approach based on an optimisation algorithm might provide moresuitable coefficients. We must also note that although the coefficients suggested in table 2
provide satisfactory predictions for several cases studied here, future research is indeed
required to examine whether these values are universal.
First, we discuss the laterally averaged velocity deficit for the Aligned Baseline (A0)case shown in figure 12. For the streamwise direction, the figure shows a sudden jump in
velocity deficit at each row, which is due to the turbine thrust force (i.e. f̄x  in (5.10)).
The model underpredicts the velocity deficit increase for the first row. This is likely due
to the lack of modelling farm-scale blockage effects. The budget analysis in § 4 showed aconsiderable flow deceleration in the vicinity of WT1 s due to farm-scale blockage effects.
Both LES data and model predictions suggest that the velocity deficit jump due to turbine
thrust forcing is significantly reduced after the first row. This is mainly due to the fact thatas shown in (5.2), the thrust force is proportional to the square of the local hub-height
velocity, which is clearly lower for subsequent turbine rows.
After the last row of turbines, the velocity deficit diminishes rather rapidly in the farmnear-wake region (e.g. for x = 50–100). The primary reason for this fast recovery in thefarm near-wake region is the large value of Reynolds stress gradient (∂u w ∂z) as shownin figure 5. The dispersive stress gradient (i.e. ∂ū ū /∂x) is large immediately behind the985 A43-24

Fast-running physics-based wake model for a wind farm

Ud = U 0 – U(a)
LES – case A0
Model – case A0
LES – case S0
Model – case S0

0.3
0.2
0.1
0
0

50

100

150

0

50

100

150

200

250

300

350

400

200

250

300

350

400

Vd = V0 – V

(b)
0.02
0
–0.02

x
Figure 12. Variation of laterally averaged (a) streamwise velocity deficit Ud and (b) spanwise velocity deficit
Vd at the turbine hub height for the Aligned Baseline (A0) and the Staggered Baseline (S0) cases. The dashedlines show the LES data and the solid lines show predictions of the developed model. Vertical dotted lines show
the locations of wind turbine rows.

ū

10
8

0.8

y 6

0.6

4
2

0.4
0

50

100

150

200

250

300

x

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 13. Contours of time-averaged streamwise velocity ū on a horizontal plane at the turbine hub height
for the Aligned Baseline (A0) case. Vertical dotted lines denote the location of wind turbine rows.

wind farm, but it decays rapidly. As discussed previously, the dispersive stress quantifies
the level of inhomogeneity in the flow. Right behind the wind farm, turbine wakes are notcompletely merged yet. This creates large velocity gradients over small length scales which
in turn lead to higher turbulence generation. Further downstream, the turbine wakes merge
and form a single farm wake. This is where the dispersive stress becomes negligible, andthe gradient of Reynolds shear stress becomes the sole mechanism for the wake recovery.
Merging of individual turbine wakes to form a single farm wake is evident in figure 13
that shows contours of the time-averaged streamwise velocity ū at the hub height. Thisis conceptually similar to the transition of a wake array to a single wake that occurs
downwind of multi-rotor turbines discussed in Bastankhah & Abkar (2019).
The developed model can capture the fast wake recovery in the farm near wake. In the
momentum equation (5.10), the streamwise velocity deficit reduction is mainly caused by
the recovery term c1 νt Ud . In the farm near wake, the streamwise velocity deficit Ud is still
fairly large. Moreover, the farm turbulent viscosity νt,f has its maximum value immediately
after the wind farm as shown in figure 10. Both of these contribute to a fast wake recovery
in the farm near-wake region. In the farm far-wake region (e.g. for x > 150), the velocity
985 A43-25

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

M. Bastankhah and others
deficit clearly decays at a slower rate. According to the budget analysis in § 4, in thisregion, the advection term is in balance with the Reynolds stress gradient, whose effect ismodelled by the recovery term c1 νt Ud in (5.10). However, both νt and Ud are much smallerin this region which leads to a slower wake recovery rate. This explains the persistence of
the wake over a very long distance. Figure 12 shows that the velocity deficit is still not
negligible even after 20 km downstream of the wind farm (i.e. x ≈ 208). It is also worth
noting that for the LES cases studied here, the Coriolis term fc Vd is an order of magnitude
smaller than other terms in the streamwise momentum equation (5.10). Moreover, while
the shear term Cx is not negligible, it is much smaller than the recovery term c1 νt Ud
for this configuration. However, this is not case for the spanwise momentum equation as
elaborated in the following.
The spanwise turbine forcing term f̄y  defined in (5.2) is zero for the LES data giventhat the yaw angle γ is zero. The Coriolis term in the spanwise momentum equation (5.10),
however, is important because it is proportional to Ud , which has a considerable value
especially within the farm. Therefore, the Coriolis term increases in the farm with anincrease of streamwise velocity deficit. According to (5.10), the Coriolis term with Ud > 0
leads to a negative Vd (i.e. positive V) in the northern hemisphere, where fc > 0. Thisis described in the literature (e.g. van der Laan & Sørensen 2017) as an anticlockwise
deflection based on a view from above. The initial anticlockwise deflection of the wake isshown in figure 12. Apart from the streamwise wake deficit, the farm-induced turbulencealso impacts the distribution of spanwise velocity deficit. In particular, the increase offarm turbulent viscosity νt,f has two important effects. It increases both the recovery term
c1 νt Vd and the veer term Cy in (5.10). Both of these effects reduce the initial anticlockwise
deflection of the wake such that after some downwind distances, the effect of the veer
term Cy becomes dominant, and the direction of the wake deflection changes to clockwise(i.e. Vd > 0). The change in the direction of wake deflection was also reported in Gadde& Stevens (2019). Ultimately, very far downstream (not shown here), Cy goes to zeroas νt,f goes to zero, and therefore the spanwise velocity deficit is completely diminishedby the recovery term. These conflicting behaviours of Coriolis and veer on the spanwise
velocity coexist as the latter is the consequence of the former. However, depending on
atmospheric conditions and simulation settings, their relative magnitude with respect
to each other could be different at different downstream positions. For instance, thesucceeding clockwise deflection observed in figure 12 can completely mask the initial
anticlockwise deflection. This may happen either in the case of a strong wind veer which
typically occurs in thermally stable boundary layers, or if the wind farm generates a highamount of turbulence. On the other hand, an anticlockwise deflection is mainly observed
if both incoming wind veer and farm-generated turbulence are relatively small, or if only
the farm near wake is of interest. This may explain why some prior works observed aclockwise deflection (e.g. Abkar et al. 2016; van der Laan & Sørensen 2017; Eriksson et al.
2019; Nouri, Vasel-Be-Hagh & Archer 2020), whereas others reported an anticlockwise
deflection (e.g. Dörenkämper et al. 2015; Allaerts & Meyers 2017; Frank et al. 2023).
Overall, the agreement of the model predictions with the LES data is satisfactory forthe streamwise velocity deficit Ud . The developed model can also successfully capture
the overall trend for the spanwise velocity deficit Vd . In agreement with the LES data, itpredicts the initial anticlockwise deflection followed by a clockwise deflection. However,
we should note that figure 12 shows some differences in the magnitude of Vd especiallyin the far wake. The value of Vd is one order of magnitude smaller than Ud for these
LES data. While predicting Vd with such small values in these cases can be challenging,the model’s prediction for Vd compared with the LES data remains within 1 % of Uh ,985 A43-26

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Fast-running physics-based wake model for a wind farmand the difference in wind direction ϕ prediction is also within 1◦ (not shown here). In
this work, we only examined neutral boundary layers where wind veer is not strong. For
thermally stable boundary layers, the cross-isobar angle θ0 in (5.28) is expected to beconsiderably larger (Peña, Gryning & Floors 2014), which in turn increases the value of
Cy . It is therefore of great interest to extend the model to thermally stratified boundary
layers in future works, where the deflection of the wake due to the wind veer is expected
to be significantly larger. The discrepancy in Vd observed between the LES data and themodel can be also partly explained by the linearisation of momentum equations. In orderto mathematically solve the equations, we replaced U with U0 = 1 in the advection term in(4.3). In the regions where the difference between U and U0 is not negligible (i.e. within
the wind farm), this assumption leads to an underestimation of the velocity deficit. The
error introduced by this assumption is expected to be more evident in spanwise velocity
predictions given their small values.
Next, we discuss the effect of wind farm layout on the farm flow distribution. The
variation of the velocity deficit for the Staggered Baseline (S0) case is shown in figure 12.
The main notable difference in the streamwise velocity deficit between the two windfarm layouts (aligned versus staggered) is the fact that, from the second row of turbines,
the jump in Ud due to turbine forcing is clearly larger for the staggered wind farm, inagreement with prior studies (e.g. Stevens et al. 2016; Stieren & Stevens 2022). This leads
to the maximum velocity deficit Ud of 0.32 for the staggered wind farm, which is 28 %larger than that for the aligned wind farm. As discussed in § 5.6, turbines within the
staggered wind farm experience a larger local velocity which according to (5.2) leads
to a larger value of turbine thrust force. This highlights the importance of the local
deficit coefficient η discussed in § 5.6 and implemented in the model (5.11). Without thiscoefficient, model predictions will be the same for both cases of A0 and S0, which isclearly unrealistic according to the LES data. The enhanced turbulence mixing caused bylarge values of velocity deficit in the staggered wind farm accelerates the wake recovery
downstream. At about x = 200, the streamwise velocity deficit becomes approximately
equal in the wake of both wind farms. The faster recovery of the wake for the staggered
wind farm is captured in the model given that the recovery term in (5.10) depends on thevalue of Ud .
The LES data also show that the spanwise velocity deficit is larger for the staggered
wind farm in comparison with the aligned wind farm. As discussed earlier, the Coriolis
term in the spanwise momentum equation (5.10) directly depends on Ud . Therefore, one
expects to observe a larger value of Vd for the staggered wind farm. Due to the initial large
anticlockwise wake deflection, the transition to the clockwise deflection occurs at a later
downstream position for the staggered wind farm. However, the wake deflection for bothwind farms eventually approaches the same value in the very far-wake region. Consistent
with the LES data, the model predicts a larger negative peak of the spanwise velocity
deficit in this case. The agreement is, however, less satisfactory for the staggered layout.
The larger streamwise velocity deficit for the S0 case may exacerbate the error introducedby the linearisation of the momentum equations discussed earlier. Nonetheless, the results
are still within a difference of 1 %Uh .
Next, we study the effect of surface roughness on the evolution of wind farm wakes.
We compare the results for the two cases of the Aligned Baseline (A0) with z0 =2 × 10−4 m/D and the Aligned Rough (AR) with z0 = 2 × 10−2 m/D. Figure 14 compares
the LES data for these two cases. Figure 14(a) shows that, as expected, the incoming
streamwise turbulence intensity I is clearly larger for the case with the higher roughness.
However, the difference in the turbulence level between the two cases becomes less
985 A43-27

(a)

1.0

0.15

U 0.9

0.10

0.8
0

20

40

60

80

100

120

140

160

180

I(%)

M. Bastankhah and others

0.05
200

LES – case A0
LES – case AR
LES – case A0
LES – case AR

I (%)
0.2

(b)
10

y

0.1

5
0

20

40

60

80

100

120

140

160

180

200

0

I (%)
0.2

(c)
10

y

0.1

5
0

20

40

60

80

100

120

140

160

180

200

0

x

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Figure 14. (a) Variation of laterally averaged streamwise turbulence intensity I and streamwise velocity U
with downwind distance x for both cases of Aligned Baseline (A0) and Aligned Rough (AR). Contours oflaterally averaged streamwise turbulence intensity I on a horizontal plane at the hub height for (b) the ARcase and (c) the A0 case.

clear within the wind farm, especially towards the end of the wind farm as shown infigure 14(a–c). In other words, these data suggest that the turbulence added by the windfarm is negatively proportional to the ambient turbulence level. This is consistent with theempirical relation of Crespo et al. (1996) for the added wake turbulence, as highlighted
in Zehtabiyan-Rezaie & Abkar (2023). In addition, it is important to note that the higher
turbulence level in the AR case promotes the wake recovery after each turbine row which
in turn increases the incoming wind speed for the next turbine. This, however, increases
the velocity deficit jump that occurs at the next turbine row according to (5.2). This is
why despite the difference in the incoming turbulence level, the two cases have fairly
similar velocity distribution within the wind farm as shown in figure 14(a). This suggeststhat the wind farm region is highly dominated by the turbulence generated by the wind
farm, and it seems to be less dependent on the incoming turbulence level. However, in thefarm wake where the turbulence added by the farm gradually diminishes, the impact of
the ambient turbulence becomes more important. Figure 14(a) shows that the wake of the
AR case recovers faster than that of the A0 case. Model predictions in comparison with
the LES data for these two cases are depicted in figure 15, which overall shows a good
agreement for the streamwise velocity deficit. The model predicts the further reduction of
the velocity deficit in the wake of the wind farm for the AR case. For the spanwise velocity
deficit, results seem to be fairly similar for the two cases. The only notable difference canbe observed in the far-wake region of the farm where the spanwise deficit for the AR caseis less than that for the A0 case. Similar to the streamwise velocity deficit, the smaller
spanwise deficit for the AR case is due to the larger value of ambient turbulence andthereby faster wake recovery.
Figures 16 and 17 show, respectively, the effect of wind farm length and turbine spacing
on the wake evolution for both the LES and the model. Figure 16 shows that the streamwise
velocity deficit is smaller for the short wind farm. This can be simply explained by the fact
985 A43-28

Fast-running physics-based wake model for a wind farm(a)
LES – case AR
Model – case AR
LES – case A0
Model – case A0

Ud = U0 – U

0.3
0.2
0.1
0
0

50

100

150

0

50

100

150

200

250

300

350

400

200

250

300

350

400

Vd = V0 – V

(b)
0.02
0
–0.02

x
Figure 15. Variation of laterally averaged (a) streamwise velocity deficit Ud and (b) spanwise velocity deficit
Vd at the turbine hub height for the Aligned Rough (AR) case. For comparison, the data for the Aligned
Baseline (A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of
the developed model. Vertical dotted lines show the locations of wind turbine rows.

(a)
LES – case AS
Model – case AS
LES – case A0
Model – case A0

0.2
0.1
0
0

50

100

150

0

50

100

150

200

250

300

350

400

200

250

300

350

400

(b)
Vd = V0 – V

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Ud = U0 – U

0.3

0.02
0
–0.02

x
Figure 16. Variation of laterally averaged (a) streamwise velocity deficit Ud and (b) spanwise velocity deficit
Vd at the turbine hub height for the Aligned Short (AS) case. For comparison, the data for the Aligned Baseline(A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of thedeveloped model. Vertical dotted lines show the locations of wind turbine rows for the AS case.

985 A43-29

M. Bastankhah and others(a)
LES – case AD
Model – case AD
LES – case A0
Model – case A0

Ud = U0 – U

0.3
0.2
0.1
0
0

50

100

150

0

50

100

150

200

250

300

350

400

200

250

300

350

400

Vd = V0 – V

(b)
0.02
0
–0.02

x
Figure 17. Variation of laterally averaged (a) streamwise velocity deficit Ud and (b) spanwise velocity deficit
Vd at the turbine hub height for the Aligned Dense (AD) case. For comparison, the data for the Aligned
Baseline (A0) case are also shown. The dashed lines show the LES data and the solid lines show predictions of
the developed model. Vertical dotted lines show the locations of wind turbine rows for the AD case.

ūh,n = 1 – ηnUd (xn)1.0
Aligned dense – case AD
0.9

Aligned baseline – case A0

0.8
0.7
0.6

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

1

2

3

4

5

6

7

8

n
Figure 18. Variation of the local hub-height velocity ūh,n as a function of turbine row number n for the
Aligned Dense (AD) and the Aligned Baseline (A0) cases based on the proposed model.

that there are fewer turbine rows in this wind farm and thereby less thrust forcing acting
on the incoming wind. In figure 17, however, the number of turbine rows are the same, andthey differ in both streamwise sx and spanwise sy inter-turbine spacing as shown in table 1.
It is interesting to note that reducing sx and sy has multifaceted effects on the laterally
averaged streamwise velocity deficit. The reduction of sy directly increases the velocity
deficit, because by reducing sy , turbine wakes occupy a larger potion of the flow field. On
the other hand, with a reduction of sx , the local incoming velocity for downwind turbines
decreases as shown in figure 18. This reduces the amount of turbine thrust force, whichis expected to decrease the total velocity deficit generated by wind turbines. However,
the turbine spacing also affects the rate of wake recovery. Decreasing both sx and sy
increases the turbulent farm velocity scale ûf , which leads to a larger νt,f and thus fasterwake recovery. Therefore, knowledge of the relative importance of each of these factors is
985 A43-30

Fast-running physics-based wake model for a wind farmneeded for each case in order to predict the overall impact of changing the turbine spacing
on the wake velocity deficit. Figure 17 shows that for this case the impact of sy changeon the velocity deficit is initially dominant as the velocity deficit within the AD farm ismuch larger than that of the A0 farm. Further downstream, however, wakes of both windfarms experience a fairly similar level of streamwise velocity deficit. The spanwise velocity
deficit is also approximately similar in both cases.

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

7. Summary

The aim of this work is to develop a new physics-based one-dimensional model to predictthe variation of laterally averaged streamwise and spanwise velocities in the wake of awind farm at the turbine hub-height level. Through a budget analysis based on the LES dataof semi-infinite wind farms, dominant terms in the momentum equations were identified.
This led to an approximate form of the momentum equations where the sum of the Coriolis
force, the divergence of the Reynolds stresses and the turbine thrust force are in balancewith the change in momentum by advection.
The linearised versions of the approximate form of the momentum equations in boththe streamwise and spanwise directions were then mathematically solved to obtain the
proposed model stated in (5.11). To derive this solution, the turbulent viscosity hypothesis
was used to model the Reynolds shear stresses. The turbulent viscosity νt was decomposed
into the ambient turbulent viscosity νt,0 and the farm turbulent viscosity νt,f (x), where thelatter changes with x. The dependency of νt,f on x was modelled using a velocity scale
proportional to the turbine forcing per unit area, and a length scale proportional to thethickness of the IBL δ. The proposed model importantly accounts for the fact that windfarms with different layouts may generate noticeably different wakes. This is mainly due
to the fact that the local incoming velocity experienced by wind turbines and thus thethrust force that is exerted by the wind turbines on the airflow depends on the farm layout.
A geometric parameter called the local deficit coefficient η was introduced to relate thelocal velocity deficit at the rotor-centre of wind turbines to the laterally averaged velocity
deficit at the same streamwise position. Moreover, the gradients of the incoming wind
shear and wind veer appeared in terms Cx and Cy of the final solution (5.11) and were
estimated using the geostrophic drag law.
The model predictions are compared with LES data for five different cases to capture theresponse of the model to various changes in farm operating conditions. The Coriolis andshear terms in the governing equation (5.10) for the streamwise direction are relatively
small. Therefore, the change in the streamwise velocity is mainly determined by how
different operating conditions influence the turbine forcing term and the wake recovery
term in (5.10). In general, a higher local incoming velocity increases the turbine forcing
term, which leads to a higher velocity deficit. The wake recovery rate on the other handis increased with an increase of either the turbulent viscosity or the velocity deficit.
Therefore, the overall impact of changing a parameter such as inter-turbine spacing or
incoming turbulence is not trivial, and these changes may lead to counteracting effects
on the streamwise velocity deficit. Overall, our results showed that the proposed model
is able to acceptably predict the variation of streamwise velocity deficit for the different
cases studied here. For the spanwise velocity deficit, turbine forces were not present in the
LES data (i.e. f̄y = 0 as γ = 0 for all turbines). However, in addition to the recovery term,
the two Coriolis and veer terms in (5.10) are of great importance. Our results showed that
the Coriolis term initially introduces an anticlockwise deflection (i.e. deflection to the left)
in the northern hemisphere. Further downstream, the veer term becomes dominant and985 A43-31

M. Bastankhah and others
introduces a clockwise deflection (i.e. deflection to the right) in the northern hemisphere.
The former is the direct effect of the Coriolis force, whereas the latter is present due to(i) an indirect effect of the Coriolis force through wind veer and (ii) the farm-generated
turbulence. The total value and direction of the wake deflection due to the Coriolis
force therefore depends on the streamwise location and more importantly on the relative
magnitude of these two counteracting terms with respect to each other.
This work serves as the first study to model the wake of a semi-infinite wind farm. While
only aligned and staggered layouts were modelled by our LESs, the developed modelis capable of modelling different layouts. Therefore, future research can implement themodel to study a wider range of layout configurations. Another interesting area of future
research is to study the effect of yaw offset on the farm wake deflection. While this is
not studied here, the effect of yaw angle is already incorporated in the model and canbe used in future works. It is also of especial interest to study how atmospheric thermal
stratification may affect the turbulent viscosity and also our estimation for the shear andveer terms. Finally, it is important to remind ourselves that by definition the model assumes
a wind farm that extends infinitely in the lateral direction. Therefore, model predictions
only resemble the flow in the wake centre of a wide finite farm where side effects are
deemed to be insignificant. More research is thus essential to quantify the impact of lateralflow entrainment and side effects. Moreover, the dimensional analysis used to simplify thevertical gradient of velocities in the recovery term of (5.10) should be scrutinised in moredetail. This simplification may imply that the cross-stream length scale associated with thevertical flow shear remains comparable to the rotor diameter. This may not be necessarily
true especially in the far wake of a wind farm.
Funding. This work was supported by Uppsala-Durham Seedcorn funding under the project entitled‘Systematic prediction of wind farm wakes: an emerging challenging in offshore wind sector’. The LESs wererun using resources provided by the Swedish National Infrastructure for Computing (SNIC).
Declaration of interests. The authors report no conflict of interest.

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Author ORCIDs.
Majid Bastankhah https://orcid.org/0000-0002-1326-027X;
Oliver R.H. Buxton https://orcid.org/0000-0002-8997-2986.
R EFERENCES
A BKAR , M. & P ORTÉ -AGEL , F. 2013 The effect of free-atmosphere stratification on boundary-layer flow andpower output from very large wind farms. Energies 6 (5), 2338–2361.
A BKAR , M. & P ORTÉ -AGEL , F. 2015 Influence of atmospheric stability on wind-turbine wakes: a large-eddy
simulation study. Phys. Fluids (1994-present) 27 (3), 035104.
A BKAR , M., S HARIFI , A. & P ORTÉ -AGEL , F. 2016 Wake flow in a wind farm during a diurnal cycle. J. Turbul.
17 (4), 420–441.
A BKAR , M., S ØRENSEN , J.N. & P ORTÉ -AGEL , F. 2018 An analytical model for the effect of vertical windveer on wind turbine wakes. Energies 11 (7), 1838.
A HSBAHS , T., N YGAARD , N.G., N EWCOMBE , A. & BADGER , M. 2020 Wind farm wakes from SAR anddoppler radar. Remote Sens. 12 (3), 462.
A LLAERTS , D. & M EYERS , J. 2017 Boundary-layer development and gravity waves in conventionally neutral
wind farms. J. Fluid Mech. 814, 95–130.
A LLAERTS , D. & M EYERS , J. 2019 Sensitivity and feedback of wind-farm-induced gravity waves. J. Fluid
Mech. 862, 990–1028.
BAI , K., K ATZ , J. & M ENEVEAU, C. 2015 Turbulent flow structure inside a canopy with complex multi-scaleelements. Boundary-Layer Meteorol. 155, 435–457.
BASTANKHAH , M. & A BKAR , M. 2019 Multirotor wind turbine wakes. Phys. Fluids 31 (8), 085106.
BASTANKHAH , M. & P ORTÉ -AGEL , F. 2014 A new analytical model for wind-turbine wakes. Renew. Energy
70 (Suppl. C), 116–123.

985 A43-32

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Fast-running physics-based wake model for a wind farm
BASTANKHAH , M. & P ORTÉ -AGEL , F. 2016 Experimental and theoretical study of wind turbine wakes inyawed conditions. J. Fluid Mech. 806, 506–541.
BASTANKHAH , M. & P ORTÉ -AGEL , F. 2017 Wind tunnel study of the wind turbine interaction with aboundary-layer flow: upwind region, turbine performance, and wake region. Phys. Fluids 29, 065105.
BASTANKHAH , M. & P ORTÉ -AGEL , F. 2019 Wind farm power optimization via yaw angle control: a windtunnel study. J. Renew. Sustain. Energy 11 (2), 023301.
BASTANKHAH , M., S HAPIRO , C.R., S HAMSODDIN , S., GAYME , D.F. & M ENEVEAU, C. 2022 A vortex
sheet based analytical model of the curled wake behind yawed wind turbines. J. Fluid Mech. 933, A2.
BASTANKHAH , M., W ELCH , B.L., M ARTÍNEZ -T OSSAS , L.A., K ING , J. & F LEMING , P. 2021 Analyticalsolution for the cumulative wake of wind turbines in wind farms. J. Fluid Mech. 911, A53.
BAY, C.J., F LEMING , P., D OEKEMEIJER , B., K ING , J., C HURCHFIELD , M. & M UDAFORT, R. 2023
Addressing deep array effects and impacts to wake steering with the cumulative-curl wake model. Wind
Energy Sci. 8 (3), 401–419.
B ELCHER , S.E., J ERRAM , N. & H UNT, J.C.R. 2003 Adjustment of a turbulent boundary layer to a canopy ofroughness elements. J. Fluid Mech. 488, 369–398.
B EMPEDELIS , N., L AIZET, S. & D ESKOS , G. 2023 Turbulent entrainment in finite-length wind farms.
J. Fluid Mech. 955, A12.
B LACKADAR , A.K. & T ENNEKES , H. 1968 Asymptotic similarity in neutral barotropic planetary boundary
layers. J. Atmos. Sci. 25 (6), 1015–1020.
B LEEG , J., P URCELL , M., RUISI , R. & T RAIGER , E. 2018 Wind farm blockage and the consequences ofneglecting its impact on energy production. Energies 11 (6), 1609.
B LONDEL , F. & C ATHELAIN , M. 2020 An alternative form of the super-gaussian wind turbine wake model.
Wind Energy Science Discuss. 2020, 1–16.
C ALAF, M., M ENEVEAU, C. & M EYERS , J. 2010 Large eddy simulation study of fully developed wind-turbine
array boundary layers. Phys. Fluids (1994-present) 22 (1), 015110.
C ALAF, M., PARLANGE , M.B. & M ENEVEAU, C. 2011 Large eddy simulation study of scalar transport infully developed wind-turbine array boundary layers. Phys. Fluids 23 (12), 126603.
C AMPAGNOLO , F., W EBER , R., S CHREIBER , J. & B OTTASSO , C.L. 2020 Wind tunnel testing of wakesteering with dynamic wind direction changes. Wind Energy Sci. 5 (4), 1273–1295.
C HAMORRO , L.P., A RNDT, R.E.A. & S OTIROPOULOS , F. 2011 Turbulent flow properties around a staggered
wind farm. Boundary-Layer Meteorol. 141 (3), 349–367.
C HAMORRO , L.P. & P ORTÉ -AGEL , F. 2009 Velocity and surface shear stress distributions behind arough-to-smooth surface transition: a simple new model. Boundary-Layer Meteorol. 130 (1), 29–41.
C HRISTIANSEN , M.B. & H ASAGER , C.B. 2005 Wake effects of large offshore wind farms identified fromsatellite SAR. Remote Sens Environ. 98 (2-3), 251–268.
C HURCHFIELD , M., L EE , S., M ORIARTY, P., M ARTINEZ , L., L EONARDI , S., V IJAYAKUMAR , G. &B RASSEUR , J. 2012a A large-eddy simulation of wind-plant aerodynamics. In 50th AIAA Aerospace
Sciences Meeting including the New Horizons Forum and Aerospace Exposition. American Institute of
Aeronautics and Astronautics.
C HURCHFIELD , M.J., L EE , S., M ICHALAKES , J. & M ORIARTY, P.J. 2012b A numerical study of the effects
of atmospheric and wake turbulence on wind turbine dynamics. J. Turbul. 13, N14.
C RESPO , A., et al. 1996 Turbulence characteristics in wind-turbine wakes. J. Wind Engng Ind. Aerodyn. 61 (1),71–85.
D ÖRENKÄMPER , M., W ITHA , B., S TEINFELD , G., H EINEMANN , D. & K ÜHN , M. 2015 The impact ofstable atmospheric boundary layers on wind-turbine wakes within offshore wind farms. J. Wind Engng Ind.
Aerodyn. 144, 146–153.
E KMAN , V.W. 1905 On the influence of the earth’s rotation on ocean-currents. Arkiv för matematik, astronomi
och fysik. 2, 11.
E RIKSSON , O., B RETON , S.P., N ILSSON , K. & I VANELL , S. 2019 Impact of wind veer and the coriolis force
for an idealized farm to farm interaction case. Appl. Sci. 9 (5), 922.
F LEMING , P., A NNONI , J., S HAH , J.J., WANG , L., A NANTHAN , S., Z HANG , Z., H UTCHINGS , K., WANG ,P., C HEN , W. & C HEN , L. 2017 Field test of wake steering at an offshore wind farm. Wind Energy Sci.
2 (1), 229–239.
F LOORS , R.R., T ROEN , I. & K ELLY, M.C. 2018 Implementation of large-scale average geostrophic wind
shear in WAsP12.1.
F RANDSEN , S. 1992 On the wind speed reduction in the center of large clusters of wind turbines. J. Wind
Engng Ind. Aerodyn. 39 (1-3), 251–265.

985 A43-33

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

M. Bastankhah and others
F RANDSEN , S.T., BARTHELMIE , R., P RYOR , S., R ATHMANN , O., L ARSEN , S., H ØJSTRUP, J. &T HØGERSEN , M. 2006 Analytical modelling of wind speed deficit in large offshore wind farms. Wind
Energy 9, 39–53.
F RANK , N.V., N EGRETTI , M.E., S OMMERIA , J., O BLIGADO , M. & C AL , R.B. 2023 Coriolis force deflects
wind plant wakes left in the northern hemisphere. arXiv:2308.01428
GADDE , S.N. & S TEVENS , R.J.A.M. 2019 Effect of coriolis force on a wind farm wake. J. Phys. Conf. Ser.
1256, 012026.
G OIT, J.P. & M EYERS , J. 2015 Optimal control of energy extraction in wind-farm boundary layers. J. Fluid
Mech. 768, 5–50.
H ANSEN , K.S., BARTHELMIE , R.J., J ENSEN , L.E. & S OMMER , A. 2012 The impact of turbulence intensity
and atmospheric stability on power deficits due to wind turbine wakes at horns rev wind farm. Wind Energy
15 (1), 183–196.
H ESS , G.D. & GARRATT, J.R. 2002 Evaluating models of the neutral, barotropic planetary boundary layerusing integral measures. Part I. Overview. Boundary-Layer Meteorol. 104, 333–358.
H OWLAND , M.F., L ELE , S.K. & DABIRI , J.O. 2019 Wind farm power optimization through wake steering.
Proc. Natl Acad. Sci. 116 (29), 14495–14500.
J ENSEN , N.O. 1983 A note on wind turbine interaction. Tech. Rep. Risø-M-2411. Risoe National Laboratory.
J ONKMAN , J., B UTTERFIELD , S., M USIAL , W. & S COTT, G. 2009 Definition of a 5-MW reference windturbine for offshore system development. Tech. Rep.. National Renewable Energy Laboratory.
K ATI Ć , I., H ØJSTRUP, J. & J ENSEN , N.O. 1986 A simple model for cluster efficiency. In Proceedings of theEuropean Wind Energy Association Conference and Exhibition, Rome, Italy, pp. 407–409.
K EANE , A., AGUIRRE , P.E.O., F ERCHLAND , H., C LIVE , P. & GALLACHER , D. 2016 An analytical model
for a full wind turbine wake. J. Phys. Conf. Ser. 753, 032039.
K ING , J., F LEMING , P., K ING , R., M ARTÍNEZ -T OSSAS , L.A., BAY, C.J., M UDAFORT, R. & S IMLEY, E.
2020 Controls-oriented model for secondary effects of wake steering. Wind Energy Sci. Discuss. 2020,
1–22.
K IRBY, A., N ISHINO , T. & D UNSTAN , T.D. 2022 Two-scale interaction of wake and blockage effects in large
wind farms. J. Fluid Mech. 953, A39.
VAN DER L AAN , M.P., BAUNGAARD , M. & K ELLY, M. 2022 Brief communication: a clarification of wakerecovery mechanisms. Wind Energy Sci. Discuss. 2022, 1–8.
VAN DER L AAN , M.P., H ANSEN , K.S., S ØRENSEN , N.N. & R ÉTHORÉ , P.E. 2015 Predicting wind farmwake interaction with RANS: an investigation of the Coriolis force. J. Phys. Conf. Ser. 625, 012026.
VAN DER L AAN , M.P., K ELLY, M., F LOORS , R. & P EÑA , A. 2020 Rossby number similarity ofan atmospheric RANS model using limited-length-scale turbulence closures extended to unstablestratification. Wind Energy Sci. 5 (1), 355–374.
VAN DER L AAN , M.P. & S ØRENSEN , N.N. 2017 Why the Coriolis force turns a wind farm wake clockwise in
the Northern Hemisphere. Wind Energy Sci. 2 (1), 285–294.
VAN DER L AAN , M.P., S ØRENSEN , N.N., R ÉTHORÉ , P.E., M ANN , J., K ELLY, M.C., T ROLDBORG , N.,
H ANSEN , K.S. & M URCIA , J.P. 2015 The k −  − fp model applied to wind farms. Wind Energy 18 (12),2065–2084.
L ANZILAO , L. & M EYERS , J. 2022 A new wake-merging method for wind-farm power prediction in thepresence of heterogeneous background velocity fields. Wind Energy 25 (2), 237–259.
L ANZILAO , L. & M EYERS , J. 2023 A parametric large-eddy simulation study of wind-farm blockage andgravity waves in conventionally neutral boundary layers. arXiv:2306.08633
L EGRIS , L., PAHUS , M.L., N ISHINO , T. & P EREZ -C AMPOS , E. 2022 Prediction and mitigation of wind farm
blockage losses considering mesoscale atmospheric response. Energies 16 (1), 386.
L ISSAMAN , P.B.S. 1979 Energy effectiveness of arbitrary arrays of wind turbines. AIAA 17, 8.
L IU, L., GADDE , S.N. & S TEVENS , R.J.A.M. 2021 Geostrophic drag law for conventionally neutral
atmospheric boundary layers revisited. Q. J. R. Meteorol. Soc. 147 (735), 847–857.
L IU, L. & S TEVENS , R.J. 2022 Vertical structure of conventionally neutral atmospheric boundary layers.
Proc. Natl Acad. Sci. 119 (22), e2119369119.
LUZZATTO -F EGIZ , P. & C OLM - CILLE , P.C. 2018 Entrainment model for fully-developed wind farms: effects
of atmospheric stability and an ideal limit for wind farm performance. Phys. Rev. Fluids 3 (9), 093802.
M ANWELL , J.F., M C G OWAN , J.G. & ROGERS , A.L. 2010 Wind Energy Explained: Theory, Design andApplication. John Wiley & Sons.
M ENEVEAU, C. 2012 The top-down model of wind farm boundary layers and its applications. J. Turbul. 13 (7),1–12.

985 A43-34

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

Fast-running physics-based wake model for a wind farm
M EYERS , J., BOTTASSO , C., DYKES , K., F LEMING , P., G EBRAAD , P., G IEBEL , G., G ÖÇMEN , T. & VAN
W INGERDEN , J.W. 2022 Wind farm flow control: prospects and challenges. Wind Energy Sci. Discuss.
2022, 1–56.
M EYERS , J. & M ENEVEAU, C. 2012 Optimal turbine spacing in fully developed wind farm boundary layers.
Wind Energy 15 (2), 305–317.
M OHAMMADI , M., BASTANKHAH , M., F LEMING , P., C HURCHFIELD , M., BOSSANYI , E., L ANDBERG , L.
& RUISI , R. 2022 Curled-skewed wakes behind yawed wind turbines subject to veered inflow. Energies
15 (23), 9135.
M OLTCHANOV, S., BOHBOT-R AVIV, Y. & S HAVIT, U. 2011 Dispersive stresses at the canopy upstream edge.
Boundary-Layer Meteorol. 139 (2), 333–351.
NARASIMHAN , G., GAYME , D.F. & M ENEVEAU, C. 2022 Effects of wind veer on a yawed wind turbine
wake in atmospheric boundary layer flow. Phys. Rev. Fluids 7 (11), 114609.
NARASIMHAN , G., GAYME , D.F. & M ENEVEAU, C. 2023 Analytical model coupling Ekman and surface
layer structure in atmospheric boundary layer flows. arXiv:2309.06650
NAVARRO D IAZ , G.P., O TERO , A.D., A SMUTH , H., S ØRENSEN , J.N. & I VANELL , S. 2022 Actuator line
model using simplified force calculation methods. Wind Energy Sci. Discuss. 2022, 1–29.
N IAYIFAR , A. & P ORTÉ -AGEL , F. 2016 Analytical modeling of wind farms: a new approach for power
prediction. Energies 9 (9), 741.
N ISHINO , T. & D UNSTAN , T.D. 2020 Two-scale momentum theory for time-dependent modelling of large
wind farms. J. Fluid Mech. 894, A2.
N OURI , R., VASEL -B E -H AGH , A. & A RCHER , C.L. 2020 The coriolis force and the direction of rotation of
the blades significantly affect the wake of wind turbines. Appl. Energy 277, 115511.
N YGAARD , N.G. & C HRISTIAN N EWCOMBE , A. 2018 Wake behind an offshore wind farm observed withdual-doppler radars. J. Phys. Conf. Ser. 1037, 072008.
N YGAARD , N.G., S TEEN , S.T., P OULSEN , L. & P EDERSEN , J.G. 2020 Modelling cluster wakes and windfarm blockage. J. Phys. Conf. Ser. 1618, 062072.
O KE , T.R. 1976 The distinction between canopy and boundary-layer urban heat islands. Atmosphere 14 (4),268–277.
P EDERSEN , J.G., SVENSSON , E., P OULSEN , L. & N YGAARD , N.G. 2022 Turbulence optimized park modelwith gaussian wake profile. J. Phys. Conf. Ser. 2265, 022063.
P EÑA , A., G RYNING , S.E. & F LOORS , R. 2014 The turning of the wind in the atmospheric boundary layer.
J. Phys. Conf. Ser. 524, 012118.
P OPE , S.B. 2000 Turbulent Flows. Cambridge University Press.
P ORTÉ -AGEL , F., BASTANKHAH , M. & S HAMSODDIN , S. 2020 Wind-turbine and wind-farm flows: a review.
Boundary-Layer Meteorol. 174 (1), 1–59.
S CHREIBER , J., BALBAA , A. & BOTTASSO , C.L. 2020 Brief communication: a double-gaussian wake model.
Wind Energy Sci. 5 (1), 237–244.
S CHUMANN , U. 1975 Subgrid scale model for finite difference simulations of turbulent flows in plane channels
and annuli. J. Comput. Phys. 18 (4), 376–404.
S COTT, R., M ARTÍNEZ -T OSSAS , L., BOSSUYT, J., H AMILTON , N. & C AL , R.B. 2023 Evolution of eddy
viscosity in the wake of a wind turbine. Wind Energy Sci. 8 (3), 449–463.
S HAPIRO , C.R., S TARKE , G.M., M ENEVEAU, C. & GAYME , D.F. 2019 A wake modeling paradigm for windfarm design and control. Energies 12 (15), 2956.
S TARKE , G.M., M ENEVEAU, C., K ING , J.R. & GAYME , D.F. 2021 The area localized coupled model for
analytical mean flow prediction in arbitrary wind farm geometries. J. Renew. Sustain. Energy 13 (3),033305.
S TEVENS , R.J. & M ENEVEAU, C. 2017 Flow structure and turbulence in wind farms. Annu. Rev. Fluid Mech.
49, 311–339.
S TEVENS , R.J.A.M., GAYME , D.F. & M ENEVEAU, C. 2015 Coupled wake boundary layer model ofwind-farms. J. Renew. Sustain. Energy 7 (2), 023115.
S TEVENS , R.J.A.M., GAYME , D.F. & M ENEVEAU, C. 2016 Effects of turbine spacing on the power output
of extended wind-farms. Wind Energy 19 (2), 359–370.
S TIEREN , A. & S TEVENS , R.J.A.M. 2022 Impact of wind farm wakes on flow structures in and arounddownstream wind farms. Flow 2, E21.
S TIPA , S., A JAY, A., A LLAERTS , D. & B RINKERHOFF , J. 2023 The multi-scale coupled model: a newframework capturing wind farm-atmosphere interaction and global blockage effects. Wind Energy Sci.
Discuss. 2023, 1–44.
S TULL , R.B. 2009 An Introduction to Boundary Layer Meteorology, vol. 13. Springer Science.
T ENNEKES , H. & LUMLEY, J.L. 1972 A First Course in Turbulence. MIT Press.

985 A43-35

M. Bastankhah and others

https://doi.org/10.1017/jfm.2024.282 Published online by Cambridge University Press

T HE MPMATH DEVELOPMENT TEAM 2023 mpmath: a Python library for arbitrary-precision floating-point
arithmetic (version 1.3.0). Available at: http://mpmath.org/
VAHIDI , D. & P ORTÉ -AGEL , F. 2022 A physics-based model for wind turbine wake expansion in theatmospheric boundary layer. J. Fluid Mech. 943, A49.
V EERS , P., et al. 2019 Grand challenges in the science of wind energy. Science 366 (6464), eaau2027.
VOUTSINAS , S., R ADOS , K. & Z ERVOS , A. 1990 On the analysis of wake effects in wind parks. Wind Engng
14, 204–219.
W HITTAKER , E.T. & WATSON , G.N. 2020 A Course of Modern Analysis. Courier Dover.
WOOD , D.H. 1982 Internal boundary layer growth following a step change in surface roughness.
Boundary-Layer Meteorol. 22 (2), 241–244.
W U, K. & P ORTÉ -AGEL , F. 2017 Flow adjustment inside and around large finite-size wind farms. Energies
10 (12), 2164.
W U, Y.T. & P ORTÉ -AGEL , F. 2012 Atmospheric turbulence effects on wind-turbine wakes: an LES study.
Energies 5 (12), 5340–5362.
X IE , S. & A RCHER , C. 2015 Self-similarity and turbulence characteristics of wind turbine wakes via
large-eddy simulation. Wind Energy 18, 1815–1838.
YANG , X., K ANG , S. & S OTIROPOULOS , F. 2012 Computational study and modeling of turbine spacing
effects in infinite aligned wind farms. Phys. Fluids 24 (11), 115107.
YOSHIZAWA , A. 1986 Statistical theory for compressible turbulent shear flows, with the application to subgrid
modeling. Phys. Fluids 29 (7), 2152–2164.
Z EHTABIYAN -R EZAIE , N. & A BKAR , M. 2023 A short note on turbulence characteristics in wind-turbine
wakes. arXiv:2302.13341
Z HANG , W., M ARKFORT, C.D. & P ORTÉ -AGEL , F. 2012 Near-wake flow structure downwind of a windturbine in a turbulent boundary layer. Exp. Fluids 52 (5), 1219–1235.
Z ONG , H. & P ORTÉ -AGEL , F. 2020 A point vortex transportation model for yawed wind turbine wakes.
J. Fluid Mech. 890, A8.

985 A43-36