A fast-running physics-based wake model for a semi-infinite wind farm

Abstract This paper presents a new generation of fast-running physics-based models to predict the wake of a semi-infinite wind farm, extending infinitely in the lateral direction but with finite 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 and spanwise directions. The developed model captures important physical phenomena such as vertical top-down transport of energy into the farm, variable wake recovery rate due to the 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 light on the counteracting effect of Coriolis force causing wake deflections in both positive and negative directions. Moreover, the impact of wind farm layout configuration on the flow 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 of various factors such as incoming turbulence, wind farm size, inter-turbine spacing and wind farm layout on the farm wake.


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 geographical areas, particularly in regions like the North Sea, where strong and consistent winds are present.Consequently, the interaction among neighboring 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 can last 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 & Meneveau 2017).Beyond technical complexities, interactions between adjacent wind farms may lead to legal and financial disputes between operators of neighboring facilities.As a result, accurate and reliable modelling of wind-farm wake effects becomes of great importance for optimizing future wind farms in increasingly competitive offshore environments.
High-fidelity numerical simulations such as large-eddy simulation (LES) are powerful tools for modeling complex turbulent wake flows, offering detailed insights into flow dynamics and wake interactions (Porté-Agel et al. 2020, and references therein).However, simulating a cluster of wind farms in congested areas such as 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 fastrunning engineering wake models striking a balance between accuracy and computational cost.Major advantages of these models are their ease of use and low computational costs, allowing for quicker assessments of various scenarios and aiding in optimization 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 modeling airflow distribution within wind farms involves predicting the wake generated by each individual turbine.A superposition method is then applied to consider the combined impact of these wake effects.These individual wake models range mainly from top-hat models (Jensen 1983;Katić et al. 1986) to Gaussian-type models (Bastankhah & Porté-Agel 2014).The Jensen top-hat model (also known as the Park model) has been extended in recent works to account for variable wake recovery rate due to turbinegenerated turbulence (Nygaard et al. 2020).Over time, Gaussian wake models have also been refined and extended in several studies to more accurately describe the near-wake region (e.g., Keane et al. 2016;Shapiro et al. 2019;Blondel & Cathelain 2020;Schreiber et al. 2020), to better capture wake expansion and its asymmetrical shape (e.g., Xie & Archer 2015;Abkar & Porté-Agel 2015;Vahidi & Porté-Agel 2022;Pedersen et al. 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 et al. 2018;Mohammadi et al. 2022;Narasimhan et al. 2022).Moreover, a variety of wake superposition methods exist, aiming to model cumulative wake effects in wind farms (e.g., Lissaman 1979;Voutsinas et al. 1990;Niayifar & Porté-Agel 2016;Zong & Porté-Agel 2020;Bastankhah et al. 2021;Lanzilao & Meyers 2022).Some of these methods are solely empirical in nature, while others have a foundation in flow physics.See Bastankhah et al. (2021) for a detailed discussion on different wake superposition methods.This simple approach has proven to be very useful in providing detailed information on the flow field within small-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 to the size of the entire wind farm or the ABL thickness.Most notably, these models fall short of capturing the crucial vertical transport of kinetic energy from higher altitude layers of the atmosphere into the wind farm/windfarm wake (Stevens & Meneveau 2017).This becomes especially problematic as the size of wind farms grows, or if we seek information about the wake of the 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 both lateral and streamwise directions, and the whole wind farm is modelled as an area with an increased aerodynamic surface roughness.Unlike single-wake modelling, this approach is able to capture the vertical transport of energy caused by turbulent fluxes, which is in balance with the energy extracted by wind turbines in infinite wind farms.The interested reader is referred to the seminal works of Frandsen (1992), Calaf et al. (2011) and other subsequent studies (Frandsen et al. 2006

Finite
-Single wake modeling -Wake superposition technique (a) 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 wind farm 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.
& Porté-Agel 2013; Stevens et al. 2016, amongst others) for more information.Despite the great advantage of these models 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 wake of the wind farm due to their core assumption that the wind farm extends infinitely in the streamwise direction.
The other group of existing models, which we classify within the broad category of the multi-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 et al. (2015) and subsequent works (e.g., Shapiro et al. 2019;Starke et al. 2021) coupled the infinite-farm approach with the single-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 et al. 2022) 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 parameterised using the classical entertainment theory (e.g., Luzzatto-Fegiz & Colm-cille 2018; Bempedelis et al. 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 in these 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 the wake 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 to perform 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 systematically model the wake of the entire farm.A schematic of the semi-infinite farm modelling in comparison with the single-turbine modelling (i.e., finite approach) as well as the infinite-farm modelling is shown in figure 1.A particular focus of this work is given to the prediction of the deflection of the farm wake.Predicting the magnitude of the farm wake deficit is important 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, the wake deflection is mainly caused by (i) meso-scale phenomena such as Coriolis force (and its by-product wind veer), and (ii) yaw misalignment.The latter has recently received a great deal of attention because of its importance in wake steering strategies (Fleming et al. 2017;Howland et al. 2019;Bastankhah & Porté-Agel 2019;Campagnolo et al. 2020, amongst others).The former, however, is mostly overlooked in prior modelling works.While the 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 have underpinned the importance of the farm-wake deflection caused by Coriolis force (e.g., Van der Laan et al. 2015;Abkar et al. 2016;Allaerts & Meyers 2017;Gadde & Stevens 2019;Eriksson et al. 2019).Interestingly, there has not been a universal agreement in the literature with regards to the direction of the wake deflection caused by the Coriolis force.van der Laan & Sørensen (2017) argued that this is due to the fact that the Coriolis force has two effects on the wake deflection.The direct effect turns the wake in the anticlockwise direction (seen from top) in the Northern hemisphere, while the indirect 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 in both streamwise and spanwise directions.The outcome is a simple one-dimensional model that predicts both laterally-averaged streamwise and spanwise velocities at the turbine hub height 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 discusses the budget analysis that is conducted to identify dominant terms in the momentum equations.The farmwake model is then developed in section 5. Results are discussed in section 6, and finally a summary is provided in section 7.

Streamwise and Spanwise Laterally-Averaged Reynolds-Averaged Navier-Stokes Equations
We start by writing the steady-state Reynolds-averaged Navier-Stokes (i.e., RANS) equation for high-Reynolds flows (i.e., negligible friction forces) using Einstein notation.For simplicity, we non-dimensionalise all variables and equations throughout this manuscript using a selection of scales based on the incoming flow and turbine characteristics.where Ω = 7.292110 −5 rad/s is the rotation rate of the earth, and  is the latitude.Note that the dimensionless Coriolis frequency   defined in (2.1) represents the ratio of Coriolis force to inertial force.This dimensionless parameter is in fact the inverse of the Rossby number  that is commonly used in geophysical studies concerning flows in oceans and atmosphere (e.g., van der Laan et al. 2020).
The dimensionless form of the RANS equation reads as (Stull 2009): where   is the velocity component in the   direction with  = 1, 2, 3 corresponding to the streamwise , spanwise  and vertical  directions, respectively.Overbar denotes time averaging, and the turbulent fluctuating velocity is  ′  =   − ū .The permutation symbol is denoted by     .Moreover, f is the time-averaged component of the turbine forces per unit volume, non-dimensionalised by U 2 ℎ /.Note that gravitational forces are neglected in (2.2) as only the streamwise and spanwise momentum directions are of interest in this study.
We separate the static pressure  in (2.2) into the one due the background driving pressure gradient   and the one due to the presence of turbines   .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 disk and its recovery to the free-stream pressure downstream.The background driving pressure gradient   can be written in terms of the geostrophic wind speed ū (Stull 2009).This simplifies (2.2) to Different forms of spatial averaging such as surface or volumetric averaging have been frequently performed in prior studies on vegetation canopies or infinite wind farms (e.g., Calaf et al. 2010;Moltchanov et al. 2011;Bai et al. 2015;Goit & Meyers 2015).In this study, however, we perform lateral averaging (averaging only along the  direction), because our aim is to determine how flow quantities at the hub-height level evolve along the streamwise direction .Mathematically, this may be defined for an arbitrary variable  as follows where ⟨⟩ indicates lateral averaging, and [−, ] is the lateral range over which the averaging is performed.The lateral fluctuation is defined as  ′′ =  − ⟨⟩ and by definition ⟨ ′′ ⟩ = 0.By performing the lateral average on (2.3), we obtain (2.5) 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 will be discussed in more details in section 4. The dispersive stress term ⟨  ′′   ′′ ⟩ that is the 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 ⟨⟩/  in (2.5) must be zero if  = 2 (i.e., / = 0).

Numerical setup: Large-eddy simulation (LES)
LESs 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 U.S. 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 the Boussinesq approximation for buoyancy effects (Churchfield et al. 2012a).
A precursor-successor approach has been utilized to develop a conventionally neutral atmospheric boundary layer 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 (U ℎ = 8 m/s) and wind direction ( = 270 • ) at the turbine hub-height level can be achieved by adjusting the magnitude and direction of the driving pressure gradient.A capping inversion with a lapse rate of 0.05 K/m is imposed at the top of the boundary layer covering the heights from 700 m to 800 m.The height at the bottom of the capping inversion, denoted by , is defined as the thickness of the ABL.The geostrophic layer above the capping inversion has a lapse rate of 0.003 K/m.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 the turbine hub-height level (Allaerts & Meyers 2017, 2019;Stipa et al. 2023;Lanzilao & Meyers 2023).The precursor simulations are run without the turbines for a period of 10 hours (36,000 seconds) to obtain a quasi-steady state.Next, the inlet conditions are recorded for a period of 9,000 seconds 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 (SGS) 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  of 126 m (Jonkman et al. 2009).The turbines are modeled as an actuator disk with no rotation and a constant thrust coefficient of 0.776.This value of   was found based on blade element momentum (BEM) simulations of a turbine rotor with U ℎ of 8 m/s and the tip-speed ratio of 7.55 (Navarro Diaz et al. 2022).The body forces are spread across the 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 to relate the average velocity on the disc with the unperturbed inflow velocity (van der Laan et al. 2015).
The same domain size is used for both precursor and successor simulations.It extends 1,000 m in the vertical direction.The first row of turbines is placed 15 rotor diameters (1,890 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 first row of turbines is relatively short in these simulations.Our aim is to maximize the use of 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., /6) long in the streamwise and lateral directions, and their height in the vertical direction grows with distance from the ground (from 2.5 m 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., /12) in the streamwise and lateral directions.This refined region starts 1,260 m (i.e., 10) 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 /24 in the streamwise and lateral directions) to capture strong velocity gradients in this 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 the wind 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 imposed for velocity and a fixed gradient for temperature.For successor simulations including wind turbines, the inlet uses the data from the precursor while a zero-gradient condition is applied at the outlet.The domain's sides, lower, and upper parts have cyclic, wall, and slip 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 in Table 1.In this table,   and   are, respectively, streamwise and spanwise inter-turbine spacing, normalised by the rotor diameter .The surface roughness normalised by  is shown by  0,0 .The number of turbine rows is denoted by , and  * is the friction velocity normalised by U ℎ .The incoming turbulence intensity at the hub height is shown by  0 , 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  for a portion of Aligned Baseline (A0) case is also shown in figure 4, where the highly-turbulent nature of the atmospheric flow and low-speed wakes are clearly visible.

Momentum Budget Analysis
In this section, the LES data for the Aligned Baseline (A0) case are employed to perform budget analysis on the laterally-averaged momentum equations (2.5).This analysis determines dominant terms in the momentum equations both within and downwind of the wind farm.This serves as a basis for the development of the physics-based model later in section 5. Note that viscous terms in the momentum equations were found to be multiple orders of 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 to steep 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 disk.4.1.Streamwise momentum equation Writing (2.5) in the streamwise direction (i.e.,  = 1) and neglecting turbine forcing gives where {, , } are respectively velocities in the streamwise, spanwise and vertical directions.
The variation of all terms in (4.1) with respect to  is illustrated in figure 5  flow quantities.As shown in figure 5 and as expected, the residual term is mostly negligible throughout the entire domain.
First, we start with the dominant streamwise advection term 1 = ⟨⟩⟨⟩/ representing the advection of streamwise momentum by the streamwise velocity.A positive value for 1 means wake recovery/flow acceleration, and vice versa.Approaching the farm, 1 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, we observe that the flow acceleration is still suppressed by farm-scale blockage effects, and 1 continues to be negative until about 3, where the maximum velocity deficit occurs (i.e., 1 = 0).The induction entrance region of the wind farm can be also illustrated by positive vertical velocity ⟨ w⟩ shown in figure 6(b).
After the second turbine row, 1 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 stress gradient 2 = ⟨ ′  ′ ⟩/ follows the profile of 1, confirming 2 acts to replenish wake momentum.In other words, peak flow acceleration (i.e., maximum 1) 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 wake recovery 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 to the common assumption that the mere presence of Reynolds stress promotes wake recovery (van der Laan et al. 2022).
Figure 5 also shows that the other advection term 2 = ⟨⟩⟨⟩/ has minimal impact on momentum transport within the domain.Moreover, although not discernible in figure 5, the Coriolis term  =   (⟨  ⟩ − ⟨⟩) is negative across the domain as expected in the northern hemisphere, but like 2, its value is negligible compared to the dominant terms.
The normal Reynolds stress gradient 1 = ⟨ ′  ′ ⟩/ is illustrative of the rate of change of turbulence level (intensity) with .1 in figure 5 indicates the turbulence level increases behind each turbine row, peaking about 3-5D downstream (i.e., where 1 = 0), before decreasing on the approach to the next row.This is in agreement with 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, 1 quickly approaches zero as the turbulence decays to its background level.
Figure 5 shows that the turbine pressure gradient  = ⟨  ⟩/ is significant within the entire farm.From actuator disk theory (Manwell et al. 2010), we know that 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 a pressure increase as wake recovery occurs.Figure 5 shows the fast recovery of pressure downwind of each turbine row indicated by positive 1.The value of 1 (i.e., rate of pressure increase) decays with  until it increases again due to the induction region of subsequent rows.It is worth noting that the variation of pressure is often neglected in wake models, 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 is even comparable to other dominant terms (e.g., 1 and 2) in the momentum equation.The term  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 in aligned wind farms (shown in figure 5) than in staggered ones (not shown here) although 1 may be still considerable within a staggered wind farm.The term 1 represents the amount of streamwise momentum transport caused by flow inhomogeneity.More precisely, it represents the rate of change in inhomogeneity with respect to the laterally-averaged flow.According to figure 5, 1 is negative within the wind farm indicating the homogeneity 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 of streamwise velocity (i.e., ū′′ ).Accordingly, the location of minimum 1 (i.e., maximum rate of approaching homogeneity) is correlated with where maximum vertical momentum transport 2 and maximum wake recovery rate 1 occur.As mixing increases so does the homogeneity, causing the reduction of the magnitude of 1 displayed in figure 5. Despite its importance within the farm, 1 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 section 6.Finally, we investigate the variation of 2 = ⟨ ′′  ′′ ⟩/.This term essentially quantifies the vertical transfer of streamwise momentum caused directly by the non-uniformity of the time-averaged windfarm flow field.Figure 5 shows that this term is clearly smaller than the other dispersive term 1.It is also interesting to note that this term is mainly positive within the wind farm.This indicates that 2 acts against wake recovery, in contrast to its turbulent counterpart 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 the wake's trajectory.Writing (2.5) in the spanwise direction (i.e.,  = 2) yields Variations of terms in (4.2) with the streamwise distance  are shown in figure 7. Due to 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 to their streamwise counterparts.
First, we start with the Coriolis term   =   (⟨⟩ − ⟨  ⟩).The dominance of this term varies across the domain.From figure 7,   is dominantly negative within the farm due to the streamwise flow deceleration, which according to (4.2) causes positive 1 = ⟨⟩⟨⟩/ and thereby positive ⟨ v⟩, as observed in figure 6(b).This deflects the wake to the left, which can be described as an anticlockwise flow rotation viewed from the top.In the wind-farm wake, figure 7 however shows that with flow acceleration the strength of the Coriolis force decays.This is where 2 = ⟨ ′  ′ ⟩/ that is the vertical turbulent entrainment of veered momentum from above becomes more dominant.Consequently, this changes the sign of the advection term 1 to negative as shown in figure 7, and therefore the far wake starts deflecting to the right (i.e., clockwise wake rotation viewed from top).In other words, the term 2 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 detail in section 6.It is also noteworthy that the magnitude of the second advection term 2, the dispersive stress terms 1 and 2 and also 1 are small, especially in the farm wake and can be neglected.

Approximate form of momentum equations
From the analysis in section 4.1, amongst others the dispersive stress (1) and pressure () 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, the summation of the four terms 1, 2, 1, and -which are challenging to model-is rather small.This is illustrated by the dashed light green color 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 1 and pressure  terms.The term 1 is negative, while  is positive; therefore, to some extent, they cancel each other out.For simplicity, we thus omit these terms from our model developed in section 5.Moreover, as discussed in section 4.2, it is apparent that in the spanwise direction the dominant terms are (i) the spanwise momentum advection by the streamwise velocity 1, (ii) the Coriolis term   and (iii) the vertical turbulent transport of veering wind 2.Therefore, the approximate forms of the momentum equations including turbine forcing can be written as (4.3) In the following section, we simplify (4.3) to develop a system of ordinary differential equations (ODEs) that can be solved mathematically for a semi-infinite wind farm.

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 {, , } aligned with the incoming wind at the turbine hub height  ℎ is adopted such that  is in the direction of the incoming wind,  represents the horizontal direction normal to , and  measures the height from the ground.The total number of wind turbine rows is denoted by .
Wind turbines in the n  ℎ row are abbreviated to WT  s, where the subscript  = {1, 2, ...,  } shows the row number labelled based on the streamwise position (i.e., ranging from  = 1 for the first row to  =  for the last row).WT  s are assumed to have the same values of thrust coefficient   , and yaw angle   , which may however be different from   , and   if  ≠ .In other words, turbines in different rows may have different operating conditions.While the lateral spacing   between turbines is assumed to be the same for all rows, the streamwise spacing between consecutive rows may be variable.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 (not done in this study).Furthermore, turbine rows may have lateral offset with respect to each other as shown in figure 8.The lateral spacing   between turbines is assumed to be constant for the whole wind farm.However, the streamwise spacing can be variable, and moreover turbine rows may be laterally shifted with respect to each other.

turbine force
The normalised aerodynamic force f exerted by wind turbines is given by where ūℎ, is the local incoming hub-height velocity for WT  s.In other words, it is the time-averaged velocity at the location of WT  s' rotor centre (i.e., (, , ) = (  ,   ,  ℎ )) in the absence of WT  s. () is the Dirac delta function, and  () is the Heaviside step function, which is defined as  () = 0 for  ⩽ 0 and  () = 1 for  > 0. We then apply the lateral averaging discussed in section 2 to obtain laterally-averaged turbine forces at the hub height  =  ℎ as follows on spatially-averaged velocity gradients, where    is the dimensionless rate of strain tensor, and   is the turbulent viscosity nondimensionalised by U ℎ .For  = 1 and  ≠ 1,   /  is an order of magnitude larger than   /  (Tennekes & Lumley 1972), so one can write (5.4) We first assess the validity of the turbulent-viscosity hypothesis using our LES data.To do so, at a given streamwise position, we examine whether −⟨ ′  ′ ⟩ is linearly proportional to ⟨ ū⟩/, according to (5.3). Figure 9 shows values of −⟨ ′  ′ ⟩ and ⟨ ū⟩/ 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 of  = 0.25 (shown in light blue) to  = 4.57 (shown in magenta).Results in figure 9 are shown for four streamwise locations, with the first two in the farm region and the other two in the wake region: (a) between WT 3 and WT 4 ( = 2.5  ), (b) between WT 7 and WT 8 ( = 6.5  ), (c) half of the farm length downstream in the wake ( = 1.5  ), where the farm length is   , (d) an entire farm length downstream in the wake ( = 2  ).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 the region 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 internal boundary layer (IBL) downwind of the wind farm due to the transition from rough to smooth terrain (see section 5.5 for more discussion on IBLs).However, further research is required to fully elucidate this phenomenon.Nevertheless, even in the farm wake, the turbulent-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   is then decomposed into two parts: one that is due to the ambient atmospheric turbulence denoted by   ,0 and one shown by   ,  corresponding to the turbulence added by the wind farm, and it varies by .In other words,   () =   ,0 +   ,  ().
(5.5) Specifying values of the ambient turbulent viscosity   ,0 and the farm turbulent viscosity   ,  is deferred to section 5.5.

Mathematical solution of velocity deficit equations
Now, we develop and solve equations for the variation of laterally-averaged velocity deficit in both streamwise   and spanwise directions   , defined as (5.6)   () =  0 −  (). (5.7) In the above equation and hereafter,  () = ⟨ ū⟩(,  =  ℎ ),  () = ⟨ v⟩(,  =  ℎ ).The subscript 0 denotes the flow in the absence of the whole wind farm.Let us recall that the coordinate system is defined based on the incoming wind direction at the hub height (see section 5.1) and all velocities in this work are non-dimensionalised by U ℎ , so  0 = 1 and  0 = 0.The latter means that   = − in this work.For generality, however, we write equations for   so the model can still be used for a different coordinate system where  0 ≠ 0. The magnitude of the local velocity deficit caused by each individual turbine can be significant especially in the near-wake region (e.g., Zhang et al. 2012;Bastankhah & Porté-Agel 2017).However, laterally-averaged values of the velocity deficit are fairly small, as will be later shown in section 6.Therefore, we linearise (4.3) by replacing ⟨ ū⟩ in the advection term with the incoming velocity  0 = 1.Moreover, using the turbulent-viscosity hypothesis discussed in section 5.3 to simplify (4.3), we obtain (5.8) In the absence of the wind farm, (5.8) is simplified to (5.9) 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 the dimensional analysis of  2   / where  1 is a constant.The term   = −  ,   2  0 / 2 is related to the rate of incoming shear, so it is called shear in (5.10).Likewise,   = −  ,   2  0 / 2 is related to the rate of incoming veer, and it is called veer in (5.10).Both   and   also depend on the amount of turbulence generated by the wind farm through   ,  .Values of   and   are determined as a function of atmospheric and farm conditions later in section 5.7, and for now they are assumed to be known.If we insert (5.2) into (5.10), the below mathematical solution for this system of ODEs can be obtained.The solution written in (5.11) is exact for a constant   ,  .For a well-defined function of   ,  () with an arbitrary distribution, this provides an approximate solution with insignificant error with respect to the exact numerical solution (not shown here) of the ODE system.
, (), (5.11) where  , () and  , () are, respectively, values of the laterally-averaged streamwise and spanwise velocity deficit caused by WT  s at , and they are given by: (5.12) In (5.12), ūℎ, is replaced with ūℎ, = 1 −     (  ), (5.13) where   (  ) =   ( =   ), and we call  the local deficit coefficient, which is defined as 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 section 5.6.To compute the ambient turbulent viscosity   ,0 and farm turbulent viscosity   ,  , (5.14) and (5.19) developed in section 5.5 are respectively used.Finally, values of   and   are estimated based on (5.30) in section 5.7.Once values of   ,   ,   ,0 and   ,  are all known, a forward marching scheme in the streamwise direction is implemented using (5.11) to find the evolution of   and   with .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.01 to 1.The solution in (5.11) is versatile as it can also be used for different distributions of   ,   and   that can potentially be developed in future studies.
5.5.Estimation of turbulent viscosity   =   ,0 +   ,  In general, turbulent viscosity can be written as a product of a turbulence velocity scale û and a turbulence length scale (i.e., mixing length) l (Tennekes & Lumley 1972).Therefore, to estimate the ambient turbulent viscosity   ,0 and the farm turbulent viscosity   ,  , one needs to specify suitable values of û and l for each term.
For the ambient turbulent viscosity   ,0 , according to Prandtl's mixing-length hypothesis (Tennekes & Lumley 1972), we assume û0 =  * and l0 =  ℎ where  ≈ 0.41 is the von-Kármań constant.Therefore, the ambient turbulent viscosity is given by   ,0 =  *  ℎ . (5.14) It is worth noting that the assumption of l0 =  ℎ is expected to be valid only in the log-law region of the ABL (Pope 2000).For our simulations the top of the rotor disk is at a location / = 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., van der Laan et al. 2020) but for simplicity we retain the simple log-law relationship for l0 .
For the farm turbulent viscosity   ,  , we use the height of the internal boundary layer (IBL) that grows above the wind farm as the turbulence length scale.Several studies have shown that the IBL grows above wind farms following the classical Elliott's  0.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   0,  = 0.28   0,  0.8

,
(5.15) where  0,  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 not account for the second IBL development.In order to use (5.15), one needs to estimate the value of  0,  for the wind farm.For an infinite wind farm, several models have been already proposed in the literature to estimate  0,  (e.g., Frandsen 1992;Calaf et al. 2010;Abkar & Porté-Agel 2013;Yang et al. 2012).The one suggested by Frandsen (1992) for an infinite wind farm with uniformly distributed wind turbines states (5.16) where    =   /(4    ),   is the normalised streamwise spacing between turbine rows, and as a reminder  0,0 is the normalised surface roughness length in the absence of the wind farm.To maintain simplicity, we use the same relationship to estimate  0,  for the semiinfinite wind farm.To compute    , we use the average streamwise inter-turbine spacing for   (i.e.,   =  −1 =1 ( +1 −   )/( − 1)), and the average value of thrust coefficient for   (i.e.,   =  =1   , /).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  (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 IBL.We therefore artificially limit the growth of the turbulence length scale l  using the below relationship, (5.17) For small values of , l  ≈ , while l  →  as the value of  → ∞.Variation of l  with  based on (5.17) is shown in figure 10.For simplicity, we here assume that the value of  is Next, we determine the turbulence velocity scale û  .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 √     0 , where  0 = 1, to estimate the turbulence velocity scale û  within the wind farm.The turbulence generation is expected to peak at some distances 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 diameter downstream (Chamorro & Porté-Agel 2009).The constant turbulence velocity scale û  caused by turbine forcing is therefore assumed to be extended from  1 to   +5 as shown in figure 10.Further downstream, û  starts to decay due to the wake recovery.The wake of a semi-infinite wind farm can be modelled as a two-dimensional wake of a canopy of roughness elements in a turbulent boundary layer.According to Belcher et al. (2003), the velocity scale defined based on the maximum velocity deficit in this type of wake flows decays with  −1 .So in summary, we model the turbulence velocity scale û  as where   = (  −  1 ) + 5. Variation of û  and   ,  are shown in figure 10, where the farm turbulent viscosity   ,  is given by (5.19) with  2 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 diameter 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   ,  for Staggered Baseline (S0) case is shown in figure 10.Fairly similar variations of turbulent viscosity have been also reported in the literature for the wake of a single wind turbine (Scott et al. 2023).

Estimation of local deficit coefficient 𝜂 𝑛
As discussed earlier, ūℎ, = 1 −     ( =   ), where   is the local deficit coefficient for WT  s.A consequence of linearising the momentum equation (5.8) is that wake effects are linearly superposed (Bastankhah et al. 2021).Therefore, one can write (5.20) where ū, (, , ) is the local velocity deficit at (, , ) caused by WT  s, and  , () is already defined in (5.12).To find ū, , one can use a Gaussian distribution exp(− 2 /2 2 ) to express the velocity deficit caused by each turbine, where  is the wake-centre velocity deficit,  is the distance from the wake centre, and  is the characteristic wake width.Therefore, for a semi-infinite wind farm, we obtain ū, where   is the maximum velocity deficit caused by each WT  at  =   ,  denotes the column number which varies from −∞ to ∞, and   is the width of WT  's wakes at  =   .Solving (5.21) gives ū, where  3 [, ] is the Jacobi theta function defined as 1 + 2 ∞ =1   2 cos(2) (Whittaker & Watson 2020).High-level programming languages (e.g., Python) may have a built-in function 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 the series.It is also noteworthy that  3 [ ± , ] =  3 [, ], so   and   in (5.22) can be the spanwise location of any turbines in WT  and WT  , respectively.
By definition  , () = ⟨ ū⟩ , (), so we can write Therefore, from (5.22) and (5.23), we conclude that ū, (5.24) The wake width   can be simplified by  * (  −   ) +  0 , where the wake expansion  * is typically around 0.02-0.04for offshore conditions and  0 is around 0.2-0.3(Bastankhah & Porté-Agel 2014).For simplicity, we assume  * = 0.025 and  0 = 0.25 to reduce (5. where  > 1 (for  = 1,  is zero), and  , () is defined in (5.12).Figure 11 shows the variation of the local deficit coefficient   with the row number  for both the Aligned Baseline (A0) and the Staggered Baseline (B0) wind farms.It is interesting to note that the value of   is always greater than one for an aligned wind farm.This is expected as in 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 and flow mixing, the value of   decays and approaches a constant value around two for this particular wind farm layout.For the staggered wind farm, however, we observe a very different behaviour.
For the second row, the value of   is almost zero.This is due to the fact that, WT 2 s are not affected by the wake of WT 1 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   .The value of   for the staggered wind farm approaches one for WT 4 and downwind rows.This suggests a fairly uniform distribution of streamwise velocity deep inside a staggered wind farm.In figure 11b, we can see how different distributions of the local velocity coefficient   lead to different distribution of local hub-height velocity ūℎ, .The local hub height velocity is clearly higher for the staggered wind farm layout which leads to more power generation.The trend shown in figure 11b is similar to the data reported in other works (e.g., Chamorro et al. 2011;Stieren & Stevens 2022).We have only discussed aligned and staggered layouts here, but an important feature of our developed model is that it is generalisable, through variable local deficit coefficient   developed in (5.26), to any conceivable wind farm layout provided the layout fulfills the requirements defined in section 5.1.As discussed earlier in (5.10),   = −  ,   2  0 / 2 and   = −  ,   2  0 / 2 .Based on (5.9),   and   can be written as (5.27) As a first-order approximation, the vertical changes in the streamwise and spanwise velocities across the ABL are proportional to  cos  0 and  sin  0 , respectively.Here,  = is the geostrophic wind speed, and the cross-isobar angle  0 is the angle between the wind direction on the ground surface and the geostrophic wind direction.One may thus write where  3 is a constant.Values of  and  0 can be obtained from the widely-used geostrophic drag law (GDL), which relates surface properties (e.g.,  0 and  * ) to geostrophic wind speed on top of the ABL (e.g., Blackadar & Tennekes 1968;Hess & Garratt 2002;van der Laan et al. 2020).For a neutrally-stratified ABL, the geostrophic drag law reads as (Liu et al. 2021) where  and  are universal empirical constants, and values of  = 1.8 and  = 4.5 used in the Wind Atlas Analysis and Application Program (WAsP) (Floors et al. 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   is positive, whereas the positive sign is for the southern hemisphere where   is negative (Liu et al. 2021).From (5.27), (5.28) and (5.29), we can therefore approximate the values of   and   as follows: (5.30)While for simplicity the conventional GDL relationship is used here to estimate   and   , recent works (e.g., Liu & Stevens 2022;Narasimhan et al. 2023) that describe the structure of Ekman boundary layer flows can be implemented in future works.

Results and discussions
In this section, we compare predictions of the model developed in section 5 with the LES data elaborated in section 3. The values of model coefficients, namely  1 ,  2 and  3 , used in this study are listed in table 2. It is worth reminding that  1 is the coefficient that is multiplied to   =   ,0 +   ,  in the final solution, while  2 is the coefficient within   ,  .An increase of either  1 or  2 increases the wake recovery rate, but  2 is the one that quantifies the importance of   ,  over   ,0 .On the other hand,  3 is the coefficient of   (i.e., shear term in the streamwise momentum equation) and   (i.e., veer term in the spanwise momentum equation).Given that   is fairly small compared to other terms in the streamwise momentum equation, the main impact of  3 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 more suitable 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 ⟩ 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 section 4 showed a considerable flow deceleration in the vicinity of WT 1 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 that as 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 farm near-wake region (e.g., for  = 50 − 100).The primary reason for this fast recovery in the farm near-wake region is the large value of Reynolds stress gradient (⟨ ′  ′ ⟩) as shown in figure 5.The dispersive stress gradient (i.e., ⟨ ′′  ′′ ⟩/) is large immediately behind the 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 not completely 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, and the 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 ū at the hub height.This is 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  1     .In the farm near wake, the streamwise velocity deficit   is still fairly large.Moreover, the farm turbulent viscosity   ,  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  > 150), the velocity deficit clearly decays at a slower rate.According to the budget analysis in section 4, in this region, the advection term is in balance with the Reynolds stress gradient, whose effect is modelled by the recovery term  1     in (5.10).However, both   and   are much smaller in 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 twenty kilometers downstream of the wind farm (i.e.,  ≈ 208).It is also worth noting that for the LES cases studied here, the Coriolis term     is an order of magnitude smaller than other terms in the streamwise momentum equation (5.10).Moreover, while the shear term   is not negligible, it is much smaller than the recovery term  1     for this configuration.However, this is not case for the spanwise momentum equation as elaborated in the following.
The spanwise turbine forcing term ⟨ f ⟩ defined in (5.2) is zero for the LES data given that the yaw angle  is zero.The Coriolis term in the spanwise momentum equation (5.10) however is important because it is proportional to   , which has a considerable value especially within the farm.Therefore, the Coriolis term increases in the farm with an increase of streamwise velocity deficit.According to (5.10), the Coriolis term with   > 0 leads to a negative   (i.e., positive ) in the northern hemisphere, where   > 0. This is described in the literature (e.g., van der Laan & Sørensen 2017) as an anticlockwise deflection based on a view from top.The initial anticlockwise deflection of the wake is shown in figure 12. Apart from the streamwise wake deficit, the farm-induced turbulence also impacts the distribution of spanwise velocity deficit.In particular, the increase of farm turbulent viscosity   ,  has two important effects.It increases both the recovery term  1     and the veer term   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   becomes dominant, and the direction of the wake deflection changes to clockwise (i.e.,   > 0).The change in the direction of wake deflection was also reported in Gadde & Stevens (2019).Ultimately, very far downstream (not shown here),   goes to zero as   ,  goes to zero, and therefore the spanwise velocity deficit is completely diminished by the recovery term.This 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, the succeeding 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 high amount 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 a clockwise deflection (e.g., Abkar et al. 2016;van der Laan & Sørensen 2017;Eriksson et al. 2019;Nouri et al. 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 for the streamwise velocity deficit   .The developed model can also successfully capture the overall trend for the spanwise velocity deficit   .In agreement with the LES data, it predicts the initial anticlockwise deflection followed by a clockwise deflection.However, we should note that figure 12 shows some differences in the magnitude of   especially in the far wake.The value of   is one order of magnitude smaller than   for these LES data.While predicting   with such small values in these cases can be challenging, the model's prediction for   compared to the LES data remains within 1% of U ℎ , and the difference in wind direction  prediction is also within 1 degree (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 be considerably larger (Peña et al. 2014), which in turn increases the value of   .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   observed between the LES data and the model can be also partly explained by the linearisation of momentum equations.In order to mathematically solve the equations, we replaced  with  0 = 1 in the advection term in (4.3).In the regions where the difference between  and  0 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 wind farm layouts (aligned vs staggered) is the fact that, from the second row of turbines, the jump in   due to turbine forcing is clearly larger for the staggered wind farm, in agreement with prior studies (e.g., Stevens et al. 2016;Stieren & Stevens 2022).This leads to the maximum velocity deficit   of 0.32 for the staggered wind farm, which is 28% larger than the one for the aligned wind farm.As discussed in section 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 section 5.6 and implemented in the model (5.11).Without this coefficient, model predictions will be the same for both cases of 0 and 0, which is clearly unrealistic according to the LES data.The enhanced turbulence mixing caused by large values of velocity deficit in the staggered wind farm accelerates the wake recovery downstream.At about  = 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 the value of   .
The LES data also shows 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   .Therefore, one expects to observe a larger value of   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 both wind 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 introduced by the linearisation of the momentum equations discussed earlier.Nonetheless, the results are still within a difference of 1%U ℎ .
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 (Case A0) with  0 = 2 × 10 −4 m/ and the Aligned Rough (Case AR) with  0 = 2 × 10 −2 m/.Figure 14 compares the LES data for these two cases.Figure 14(a) shows that, as expected, the incoming streamwise turbulence intensity ⟨⟩ is clearly larger for the case with the higher roughness.However, the difference in the turbulence level between the two cases become less clear within the wind farm, especially towards the end of the wind farm as shown in figure 14(a-c).In other words, these data suggest that the turbulence added by the wind farm is negatively proportional to the ambient turbulence level.This is consistent with the empirical 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 14(a).This suggests that the wind farm region is highly dominated by the turbulence generated by the wind farm, and it seems to be less dependant on the incoming turbulence level.However, in the farm 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 the one 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 can be observed in the far-wake region of the farm where the spanwise deficit for the AR case is less than the one for 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 and thereby 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 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, and they differ in both streamwise   and spanwise   inter-turbine spacing as shown in table 1.It is interesting to note that reducing   and   has multifaceted effects on the laterally-averaged streamwise velocity deficit.The reduction of   directly increases the velocity deficit, because by reducing   , turbine wakes occupy a larger potion of the flow field.On the other hand, with a reduction of   , the local incoming velocity for downwind turbines decreases as shown in figure 18.This reduces the amount of turbine thrust force, which is expected to decrease the total velocity deficit generated by wind turbines.However, the turbine spacing also affects the rate of wake recovery.Decreasing both   and   increases the turbulent farm velocity scale û  , which leads to a larger   ,  and thus faster wake recovery.Therefore, knowledge on the relative importance of each of these factors is needed 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   change on the velocity deficit is initially dominant as the velocity deficit within the AD farm is much larger than the one of the A0 farm.Further downstream, however, wakes of both wind farms experience a fairly similar level of streamwise velocity deficit.The spanwise velocity deficit is also approximately similar in both cases.

Summary
The aim of this work is to develop a new physics-based one-dimensional model to predict the variation of laterally-averaged streamwise and spanwise velocities in the wake of a wind farm at the turbine hub-height level.Through a budget analysis based on the LES data of 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 balance with the change in momentum by advection.
The linearised versions of the approximate form of the momentum equations in both the streamwise and spanwise directions were then mathematically solved to obtain the 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   was decomposed into the ambient turbulent viscosity   ,0 and the farm turbulent viscosity   ,  (), where the latter changes with .The dependency of   ,  to  was modelled using a velocity scale proportional to the turbine forcing per unit area, and a length scale proportional to the thickness of the internal boundary layer .The proposed model importantly accounts for the fact that wind farms 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 the thrust 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 the local 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   and   of the final solution (5.11) and were estimated using the Geostrophic drag law.
The model predictions are compared to LES data for five different cases to capture the response of the model to various changes in farm operating conditions.The Coriolis and shear 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 hand is 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 = 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 and 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 model is capable of modelling different layouts.Therefore, future research can implement the model 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 can be used in future works.It is also of special interest to study how atmospheric thermal stratification may affect the turbulent viscosity and also our estimation for the shear and veer 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 lateral flow entrainment and side effects.Moreover, the dimensional analysis used to simplify the vertical gradient of velocities in the recovery term of (5.10) should be scrutinised in more detail.This simplification may imply that the cross-stream length-scale associated with the vertical flow shear remains comparable to the rotor diameter.This may not be necessarily true especially in the far wake of a wind farm.

-
Infinite width, finite length -Predict streamwise and spanwise deficit

Figure 2 :
Figure 2: Schematic of the computational domain for the A0 case.Turbines are shown by black circles, and the rotor diameter is denoted by .
turbine spacing, normalised by the rotor diameter .The surface roughness normalised by  is shown by  0,0 .The number of turbine rows is denoted by , and  * is the friction velocity normalised by U ℎ .The incoming turbulence intensity at the hub height is shown by  0 , 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).

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

Figure 5 :
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-dimensionalised using a selection of U ℎ and .

Figure 6 :
Figure 6: Variation of laterally-averaged (a) streamwise ⟨ ū⟩, (b) spanwise ⟨ v⟩ and vertical ⟨ w⟩ velocities with .All variables are non-dimensionalised using a selection of U ℎ and .The locations of turbine rows are denoted by vertical dotted lines.

Figure 7 :
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-dimensionalised using a selection of U ℎ and .

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

Figure 9 :
Figure 9: Distribution of −⟨ ′  ′ ⟩ against ⟨ ū⟩/ at different heights for four different streamwise positions in the A0 case, where (a-b) are in the farm region, and (c-d) are in the wake region, and   = 7  is the farm length.

Figure 10 :
Figure 10: Variation of velocity scale û  , length scale l  and farm turbulent viscosity   , based on the modelling approached elaborated in section 5.5.Variation of   ,  for the LES data (case S0) is also shown.In the figure, û  is normalised by √    ,   ,  is normalised by   ,  ( =  1 +   ), and l  is normalised by .Vertical dotted lines show the location of turbine rows, and the vertical dashed line shows  1 +   where the turbulent viscosity   ,  is maximum.

Figure 11 :
Figure 11: Variation of (a) local deficit coefficient   and (b) local hub-height velocity ūℎ, with the row number  for the Staggered Baseline (S0) and the Aligned baseline (A0) cases based on (5.26).

Figure 12 :
Figure 12: Variation of laterally-averaged streamwise velocity deficit   and spanwise velocity deficit   at the turbine hub height for the Aligned Baseline (A0) and the Staggered Baseline (S0) case.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.

Figure 13 :
Figure 13: Contours of time-averaged streamwise velocity ū on a horizontal plane at the turbine hub height for Aligned Baseline (A0) case.Vertical dotted lines denote the location of wind turbine rows.
Figure 14: (a) Variation of laterally-averaged streamwise turbulence intensity ⟨⟩ and streamwise velocity  with downwind distance  for both cases of Aligned Baseline (A0) and Aligned Rough (AR).Contours of laterally-averaged streamwise turbulence intensity ⟨⟩ on a horizontal plane at the hub height for (b) the AR case and (c) for the A0 case.

Figure 15 :
Figure 15: Variation of the laterally-averaged streamwise velocity deficit   and the spanwise velocity deficit   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.

Figure 16 :
Figure 16: Variation of laterally-averaged streamwise velocity deficit   and spanwise velocity deficit   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 the developed model.Vertical dotted lines show the locations of wind turbine rows for the AS case.

Figure 17 :
Figure 17: Variation of laterally-averaged streamwise velocity deficit   and spanwise velocity deficit   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.

Figure 18 :
Figure 18: Variation of the local hub-height velocity ūℎ, as a function of turbine row number  for the Aligned Dense (AD) and the Aligned Baseline (A0) cases based on the proposed model.
2 ℎ , where  is the air density.The dimensionless Coriolis frequency   is defined as All spatial dimensions are normalised by the turbine rotor diameter .All velocities   are normalised by the incoming velocity U ℎ at the turbine hub height  ℎ .Static pressure  is normalised by Focus on Fluids articles must not exceed this page length U

Table 2 :
Empirical model coefficients used in this study.5.7.Estimation of shear term   and veer term