4.1. Streamwise momentum equation

Writing (2.5) in the streamwise direction (i.e. $i=1$) and neglecting turbine forcing gives

(4.1)\begin{align} \underbrace{\langle\bar{u}\rangle\frac{\partial\langle\bar{u}\rangle}{\partial x}}_{{A1x}} + \underbrace{\langle\bar{w}\rangle \frac{\partial\langle\bar{u}\rangle}{\partial z}}_{{A2x}} &={-}\underbrace{f_c(\langle\overline{v_g}\rangle - \langle\bar{v}\rangle)\vphantom{\frac{\partial}{\partial}}}_{{Cx}} - \underbrace{\frac{\partial\langle\bar{u}''\bar{u}''\rangle}{\partial x}}_{{D1x}}\nonumber\\ &\quad - \underbrace{\frac{\partial\langle\bar{u}''\bar{w}''\rangle}{\partial z}}_{{D2x}} - \underbrace{\frac{\partial\langle\overline{u'u'}\rangle}{\partial x}}_{{R1x}} - \underbrace{\frac{\partial\langle\overline{u'w'}\rangle}{\partial z}}_{{R2x}} - \underbrace{\frac{\partial \langle \bar{p}_t\rangle}{\partial x}}_{{Px}}, \end{align}

where $\{u,v,w\}$ are respectively velocities in the streamwise, spanwise and vertical directions. 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 is observed in the flow quantities. As shown in figure 5 and as expected, the residual term is mostly negligible throughout the entire domain.

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 $\mathcal {U}_h$ and $D$.

First, we start with the dominant streamwise advection term $A1x=\langle \bar {u}\rangle {\partial \langle \bar {u}\rangle }/{\partial 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. Reference Bleeg, Purcell, Ruisi and Traiger2018). Behind the first row of turbines, we observe that the flow acceleration is still suppressed by farm-scale blockage effects, and $A1x$ 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 by positive vertical velocity $\langle \bar {w}\rangle$ shown in figure 6(b).

Figure 6. Variation of laterally averaged (a) streamwise $\langle \bar {u}\rangle$ and (b) spanwise $\langle \bar {v}\rangle$ and vertical $\langle \bar {w}\rangle$ velocities with $x$. All variables are non-dimensionalised using a selection of $\mathcal {U}_h$ and $D$. The locations of turbine rows are denoted 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 stress gradient $R2x={\partial \langle \overline {u'w'}\rangle }/{\partial 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 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, Baungaard & Kelly Reference van der Laan, Baungaard and Kelly2022).

Figure 5 also shows that the other advection term $A2x=\langle \bar {w}\rangle {\partial \langle \bar {u}\rangle }/{\partial z}$ has minimal impact on momentum transport within the domain. Moreover, although not discernible in figure 5, the Coriolis term $Cx=f_c(\langle \overline {v_g}\rangle - \langle \bar {v}\rangle )$ is negative across the domain as expected in the northern hemisphere, but like $A2x$, its value is negligible compared with the dominant terms.

The normal Reynolds stress gradient $R1x={\partial \langle \overline {u'u'}\rangle }/\partial x$ is illustrative of the rate of change 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. where $R1x=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 Reference Wu and Porté-Agel2012). 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={\partial \langle \bar {p}_t\rangle }/{\partial x}$ is significant within the entire farm. From actuator disc theory (Manwell, McGowan & Rogers Reference Manwell, McGowan and Rogers2010), 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 $P1x$. The value of $P1x$ (i.e. rate of pressure increase) decays with $x$ 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. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021) 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. $A1x$ and $R2x$) in the momentum equation. 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. Reference Moltchanov, Bohbot-Raviv and Shavit2011). 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. Reference Moltchanov, Bohbot-Raviv and Shavit2011). 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 $D1x$ may be still considerable within a staggered wind farm. The term $D1x$ 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, $D1x$ 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. $\bar {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={\partial \langle \bar {u}''\bar {w}''\rangle }/{\partial 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 term is clearly smaller than the other dispersive term $D1x$. It is also interesting to note that this term is mainly positive within the wind farm. This indicates that $D2x$ acts against wake recovery, 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 the wake's trajectory. Writing (2.5) in the spanwise direction (i.e. $i=2$) yields

(4.2)\begin{align} \underbrace{\langle\bar{u}\rangle\frac{\partial\langle\bar{v}\rangle}{\partial x}}_{{A1y}} + \underbrace{\langle\bar{w}\rangle\frac{\partial\langle\bar{v}\rangle}{\partial z}}_{{A2y}} &={-}\underbrace{f_c(\langle\bar{u}\rangle-\langle\overline{u_g}\rangle) \vphantom{\frac{\partial}{\partial}}}_{{Cy}} - \underbrace{\frac{\partial \langle\bar{v}''\bar{u}''\rangle}{\partial x}}_{{D1y}} - \underbrace{\frac{\partial\langle\bar{v}''\bar{w}''\rangle}{\partial z}}_{{D2y}}\nonumber\\ &\quad -\underbrace{\frac{\partial\langle\overline{v'u'}\rangle}{\partial x}}_{{R1y}} - \underbrace{\frac{\partial\langle\overline{v'w'}\rangle}{\partial z}}_{{R2y}}. \end{align}

Variations of terms in (4.2) with the streamwise distance $x$ 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 with their streamwise counterparts.

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 $\mathcal {U}_h$ and $D$.

First, we start with the Coriolis term $Cy=f_c(\langle \bar {u}\rangle -\langle \overline {u_g}\rangle )$. The dominance of this term varies across the domain. From figure 7, $Cy$ is dominantly negative within the farm due to the streamwise flow deceleration, which according to (4.2) causes positive $A1y=\langle \bar {u}\rangle {\partial \langle \bar {v}\rangle }/{\partial x}$ and thereby positive $\langle \bar {v} \rangle$, 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 $R2y={\partial \langle \overline {v'w'}\rangle }/{\partial z}$ that is the vertical turbulent entrainment of veered momentum from above becomes more dominant. Consequently, this changes the sign of the advection term $A1y$ 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 $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 detail in § 6. It is also noteworthy that the magnitudes of the second advection term $A2y$, the dispersive stress terms $D1y$ and $D2y$ and also $R1y$ are small, especially in the farm wake and 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, the summation of the four terms $D1x$, $D2x$, $R1x$ and $Px$ – which are challenging to model – is rather 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 developed in § 5. Moreover, as discussed in § 4.2, it is apparent that in the spanwise direction the dominant 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 can be written as

(4.3a)$$\begin{gather} \langle \bar{u}\rangle \frac{\partial \langle \bar{u}\rangle}{\partial x} \approx{-}f_c (\langle \overline{v_{\mathit{g}}}\rangle -\langle \bar{v}\rangle)\vphantom{\frac{\rangle\partial}{\partial}} - \frac{\partial \langle \overline{u'w'}\rangle}{\partial z} -\langle\, \bar{f}_{x} \rangle, \end{gather}$$

(4.3b)$$\begin{gather}\langle \bar{u}\rangle \frac{\partial \langle \bar{v}\rangle}{\partial x} \approx{+}f_c (\langle \overline{u_{\mathit{g}}}\rangle -\langle \bar{u}\rangle)\vphantom{\frac{\rangle\partial}{\partial}} - \frac{\partial \langle \overline{v'w'}\rangle}{\partial z} -\langle\, \bar{f}_{y} \rangle. \end{gather}$$

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.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 $z_h$ is adopted such that $x$ is in the direction of the incoming wind, $y$ represents the horizontal direction normal to $x$ and $z$ measures the height from the ground. The total number of wind turbine rows is denoted by $N$. Wind turbines in the $n$th row are abbreviated to WT$_n$s, where the subscript $n=\{1,2,\ldots,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 WT$_n$s are assumed to have the same values of thrust coefficient $C_{T,n}$ and yaw angle $\gamma _n$, which may, however, be different from $C_{T,m}$ and $\gamma _m$ if $n\neq m$. In other words, turbines in different rows may have different operating conditions. While the lateral spacing $s_y$ 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 position of WT$_n$s is $y_n+ks_y$, where $k=-\infty \dots \infty$.

Figure 8. Schematic of a semi-infinite wind farm. Turbines in row $n$ are denoted by WT$_n$. The lateral spacing $s_y$ 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.

5.2. Turbine force

The normalised aerodynamic force $\bar {f}$ exerted by wind turbines is given by

(5.1a)$$\begin{gather} \bar{f}_x ={-}\sum_{n=1}^N\sum_{k={-}\infty}^\infty\frac{1}{2} C_{T,n}\bar{u}_{h,n}^2\cos(\gamma_n)\delta(x-x_n)H(0.5^2-[(y-y_n-ks_y)^2+(z-z_h)^2]), \end{gather}$$

(5.1b)$$\begin{gather}\bar{f}_y={-}\sum_{n=1}^N\sum_{k={-}\infty}^\infty\frac{1}{2} C_{T,n}\bar{u}_{h,n}^2\sin(\gamma_n)\delta(x-x_n)H(0.5^2-[(y-y_n-ks_y)^2+(z-z_h)^2]), \end{gather}$$

where $\bar {u}_{h,n}$ is the local incoming hub-height velocity for WT$_n$s. In other words, it is the time-averaged velocity at the location of WT$_n$s’ rotor centre (i.e. $(x,y,z)=(x_n,y_n,z_h)$) in the absence of WT$_n$s. Here $\delta (x)$ is the Dirac delta function and $H(x)$ is the Heaviside step function, which is defined as $H(x)=0$ for $x\leq 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 height $z=z_h$ as follows:

(5.2a)$$\begin{gather} \langle\,\bar{f}_x\rangle ={-}\frac{1}{2s_y}\sum_{n=1}^N C_{T,n}\bar{u}_{h,n}^2\cos(\gamma_n)\delta(x-x_n), \end{gather}$$

(5.2b)$$\begin{gather}\langle\,\bar{f}_y\rangle ={-}\frac{1}{2s_y}\sum_{n=1}^N C_{T,n}\bar{u}_{h,n}^2\sin(\gamma_n)\delta(x-x_n). \end{gather}$$

5.3. Simplifying Reynolds shear stress terms in (4.3)

The Boussinesq eddy-viscosity hypothesis has been used in previous studies (Belcher, Jerram & Hunt Reference Belcher, Jerram and Hunt2003; Luzzatto-Fegiz & Colm-cille Reference Luzzatto-Fegiz and Colm-cille2018) to model spatially averaged Reynolds stresses based on spatially averaged velocity gradients:

(5.3)\begin{equation} \langle \overline{u_i'u_j'}\rangle ={-}2\nu_t \left[ \underbrace{\frac{1}{2}\left(\frac{\partial \langle \overline{u_i}\rangle}{\partial x_j} + \frac{\partial \langle \overline{u_j}\rangle}{\partial x_i}\right)}_{{S_{ij}}}\right]\quad (\textrm{for}\ i\neq j), \end{equation}

where $S_{ij}$ is the dimensionless rate of strain tensor and $\nu _t$ is the turbulent viscosity non-dimensionalised by $\mathcal {U}_hD$. For $i=1$ and $j\neq 1$, $\partial u_i/\partial x_j$ is an order of magnitude larger than $\partial u_j/\partial x_i$ (Tennekes & Lumley Reference Tennekes and Lumley1972), so one can write

(5.4a,b)\begin{equation} \frac{\partial \langle \overline{u'w'}\rangle}{ \partial z} \approx{-}\nu_t \frac{\partial^2 \langle \bar{u}\rangle}{\partial z^2},\quad \frac{\partial \langle \overline{v'w'}\rangle}{\partial z} \approx{-}\nu_t \frac{\partial^2 \langle \bar{v}\rangle}{\partial z^2}. \end{equation}

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 $-\langle \overline {u'w'}\rangle$ is linearly proportional to $\partial \langle \bar {u}\rangle /\partial z$, according to (5.3). Figure 9 shows values of $-\langle \overline {u'w'}\rangle$ and $\partial \langle \bar {u}\rangle /\partial 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 of $z=0.25$ (shown in light blue) to $z=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: between WT$_3$ and WT$_4$ ($x=2.5s_x$) (figure 9a), between WT$_7$ and WT$_8$ ($x=6.5s_x$) (figure 9b), half of the farm length downstream in the wake ($x=1.5x_N$), where the farm length is $x_N$ (figure 9c) and an entire farm length downstream in the wake ($x=2x_N$) (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 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 § 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.

Figure 9. Distribution of $-\langle \overline {u'w'}\rangle$ against $\partial \langle \bar {u}\rangle /\partial z$ 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 $x_N=7s_x$ is the farm length.

The turbulent viscosity $\nu _t$ is then decomposed into two parts: one that is due to the ambient atmospheric turbulence denoted by $\nu _{t,0}$ and one shown by $\nu _{t,f}$ corresponding to the turbulence added by the wind farm, and it varies by $x$. In other words,

(5.5)\begin{equation} \nu_t(x)=\nu_{t,0}+\nu_{t,f}(x). \end{equation}

Specifying values of the ambient turbulent viscosity $\nu _{t,0}$ and the farm turbulent viscosity $\nu _{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 $U_d$ and spanwise $V_d$ directions, defined as

(5.6)$$\begin{gather} U_d(x)=U_0- U(x), \end{gather}$$

(5.7)$$\begin{gather}V_d(x)=V_0-V(x). \end{gather}$$

In the above equation and hereafter, $U(x)=\langle \bar {u}\rangle (x,z=z_h)$, $V(x)=\langle \bar {v}\rangle (x,z=z_h)$. 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 § 5.1) and all velocities in this work are non-dimensionalised by $\mathcal {U}_h$, so $U_0=1$ and $V_0=0$. The latter means that $V_d=-V$ in this work. For generality, however, we write equations for $V_d$ so the model can still be used for a different coordinate system where $V_0\neq 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, Markfort & Porté-Agel Reference Zhang, Markfort and Porté-Agel2012; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2017). However, laterally averaged values of the velocity deficit are fairly small, as is later shown in § 6. Therefore, we linearise (4.3) by replacing $\langle \bar {u}\rangle$ in the advection term with the incoming velocity $U_0=1$. Moreover, using the turbulent-viscosity hypothesis discussed in § 5.3 to simplify (4.3), we obtain

(5.8a)$$\begin{gather} \frac{\partial U}{\partial x} \approx{-}f_c (V_{\mathit{g}} -V) +\nu_t\frac{\partial^2 U}{\partial z^2} -\langle\, \bar{f}_{x} \rangle, \end{gather}$$

(5.8b)$$\begin{gather}\frac{\partial V}{\partial x} \approx{+}f_c (U_{\mathit{g}} -U) +\nu_t \frac{\partial^2 V}{\partial z^2} -\langle\, \bar{f}_{y} \rangle. \end{gather}$$

In the absence of the wind farm, (5.8) is simplified to

(5.9a)$$\begin{gather} 0 \approx{-}f_c (V_{\mathit{g}} -V_0)\vphantom{\frac{\rangle\partial}{\partial}} +\nu_{t,0} \frac{\partial^2 U_0}{\partial z^2}, \end{gather}$$

(5.9b)$$\begin{gather}0 \approx{+}f_c (U_{\mathit{g}} -U_0) +\nu_{t,0} \frac{\partial^2 V_0}{\partial z^2}. \end{gather}$$

It is worth noting that (5.9) was solved in Ekman (Reference Ekman1905) 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 $\partial ^2 U_d/\partial z^2\propto U_d/l^2$, where $l\approx 1$ (i.e. length scale comparable to rotor diameter), we obtain

(5.10a)$$\begin{gather} \frac{\textrm{d} U_d}{\textrm{d} x} \approx{+} \underbrace{f_cV_d\vphantom{\frac{a}{b}}}_{Coriolis} -\underbrace{c_1 \nu_t U_d \vphantom{\frac{a}{b}}}_{recovery} +\underbrace{\langle\, \bar{f}_{x}\rangle \vphantom{\frac{a}{b}}}_{turbine\ forcing}+ \underbrace{C_{x}\vphantom{\frac{a}{b}}}_{shear}, \end{gather}$$

(5.10b)$$\begin{gather}\frac{\textrm{d} V_d}{\textrm{d} x} \approx{-} \underbrace{f_c U_d\vphantom{\frac{a}{b}}}_{Coriolis} -\underbrace{c_1 \nu_t V_d \vphantom{\frac{V_d}{Ro}}}_{recovery} +\underbrace{\langle\, \bar{f}_{y}\rangle \vphantom{\frac{V_d}{Ro}}}_{turbine\ forcing}+ \underbrace{C_{y}\vphantom{\frac{V_d}{Ro}}}_{veer}, \end{gather}$$

where $c_1$ is a constant. The term $C_{x}=-\nu _{t,f} \partial ^2 U_0/\partial z^2$ is related to the rate of incoming shear, so it is called shear in (5.10). Likewise, $C_{y}=-\nu _{t,f} \partial ^2 V_0/\partial z^2$ is related to the rate of incoming veer, and it is called veer in (5.10). Both $C_x$ and $C_y$ also depend on the amount of turbulence generated by the wind farm through $\nu _{t,f}$. Values of $C_x$ and $C_y$ are determined as a 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 system of ordinary differential equations can be obtained. The solution written in (5.11) is exact for a constant $\nu _{t,f}$. For a well-defined function of $\nu _{t,f}(x)$ with an arbitrary distribution, this provides an approximate solution with insignificant error with respect to the exact numerical solution (not shown here) of the ordinary differential equation system:

(5.11a)$$\begin{gather} U_d(x)\approx \frac{C_x}{c_1\nu_t }(1-\textrm{e}^{{-}c_1\int_{0}^{x} \nu_t\,\textrm{d} x})H(x)+\sum_{n=1}^N U_{d,n}(x), \end{gather}$$

(5.11b)$$\begin{gather}V_d(x)\approx \frac{C_y}{c_1\nu_t }(1-\textrm{e}^{{-}c_1\int_{0}^{x} \nu_t\,\textrm{d} x})H(x)+\sum_{n=1}^N V_{d,n}(x), \end{gather}$$

where $U_{d,n}(x)$ and $V_{d,n}(x)$ are, respectively, values of the laterally averaged streamwise and spanwise velocity deficit caused by WT$_n$s at $x$, and they are given by

(5.12a)$$\begin{gather} U_{d,n}(x)=\frac{C_{T,n}}{2s_y}(1-\eta_n U_d(x_n))^2\,\textrm{e}^{{-}c_1\int_{x_n}^{x} \nu_t\,\textrm{d} x}\cos(\gamma_n-f_c(x-x_n))H(x-x_n), \end{gather}$$

(5.12b)$$\begin{gather}V_{d,n}(x)=\frac{C_{T,n}}{2s_y}(1-\eta_n U_d(x_n))^2\,\textrm{e}^{{-}c_1\int_{x_n}^{x} \nu_t\,\textrm{d} x}\sin(\gamma_n-f_c(x-x_n))H(x-x_n). \end{gather}$$

In (5.12), $\bar {u}_{h,n}$ is replaced with

(5.13)\begin{equation} \bar{u}_{h,n}=1-\eta_n U_d(x_n), \end{equation}

where $U_d(x_n)=U_d(x=x_n)$, and we call $\eta$ 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 § 5.6. To compute the ambient turbulent viscosity $\nu _{t,0}$ and farm turbulent viscosity $\nu _{t,f}$, (5.14) and (5.19) developed in § 5.5 are respectively used. Finally, values of $C_x$ and $C_y$ are estimated based on (5.30) in § 5.7. Once values of $C_x$, $C_y$, $\nu _{t,0}$ and $\nu _{t,f}$ are all known, a forward marching scheme in the streamwise direction is implemented using (5.11) to find the evolution of $U_d$ and $V_d$ with $x$. 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 $\nu _{t}$, $C_x$ and $C_y$ that can potentially be developed in future studies.

5.5. Estimation of turbulent viscosity $\nu _t=\nu _{t,0}+\nu _{t,f}$

In general, turbulent viscosity can be written as a product of a turbulence velocity scale $\hat {u}$ and a turbulence length scale (i.e. mixing length) $\hat {l}$ (Tennekes & Lumley Reference Tennekes and Lumley1972). Therefore, to estimate the ambient turbulent viscosity $\nu _{t,0}$ and the farm turbulent viscosity $\nu _{t,f}$, one needs to specify suitable values of $\hat {u}$ and $\hat {l}$ for each term.

For the ambient turbulent viscosity $\nu _{t,0}$, according to Prandtl's mixing-length hypothesis (Tennekes & Lumley Reference Tennekes and Lumley1972), we assume $\hat {u}_0=u_*$ and $\hat {l}_0=\kappa z_h$, where $\kappa \approx 0.41$ is the von Kármań constant. Therefore, the ambient turbulent viscosity is given by

(5.14)\begin{equation} \nu_{t,0}=\kappa u_*z_h. \end{equation}

It is worth noting that the assumption of $\hat {l}_0=\kappa z_h$ is expected to be valid only in the log-law region of the ABL (Pope Reference Pope2000). 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. van der Laan et al. Reference van der Laan, Kelly, Floors and Peña2020) but for simplicity we retain the simple log-law relationship for $\hat {l}_0$.

For the farm turbulent viscosity $\nu _{t,f}$, we use the height of the 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 $x^{0.8}$ power law (e.g. Allaerts & Meyers Reference Allaerts and Meyers2017; Wu & Porté-Agel Reference Wu and Porté-Agel2017). According to Wood (Reference Wood1982), the thickness of the IBL, $\delta$, due to a smooth to rough transition is given by

(5.15)\begin{equation} \frac{\delta}{z_{0,f}}=0.28\left(\frac{x}{z_{0,f}}\right)^{0.8}, \end{equation}

where $z_{0,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 Reference Oke1976). 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 $z_{0,f}$ for the wind farm. For an infinite wind farm, several models have been already proposed in the literature to estimate $z_{0,f}$ (e.g. Frandsen Reference Frandsen1992; Calaf et al. Reference Calaf, Meneveau and Meyers2010; Yang et al. Reference Yang, Kang and Sotiropoulos2012; Abkar & Porté-Agel Reference Abkar and Porté-Agel2013). The one suggested by Frandsen (Reference Frandsen1992) for an infinite wind farm with uniformly distributed wind turbines states

(5.16)\begin{equation} z_{0,f}=z_h\exp\left(-\dfrac{\kappa}{\sqrt{\dfrac{1}{2}c_{ft} +\left[\dfrac{\kappa}{\ln(z_h/z_{0,0})}\right]^2}}\right), \end{equation}

where $c_{ft}={\rm \pi} C_T/(4s_xs_y)$, $s_x$ is the normalised streamwise spacing between turbine rows and as a reminder $z_{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 $z_{0,f}$ for the semi-infinite wind farm. To compute $c_{ft}$, we use the average streamwise inter-turbine spacing for $s_x$ (i.e. $s_x=\sum _{n=1}^{N-1}(x_{n+1}-x_{n})/(N-1)$), and the average value of thrust coefficient for $C_T$ (i.e. $C_T=\sum _{n=1}^{N}C_{T,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 $\delta$ is smaller than the ABL thickness $H$ (Wood Reference Wood1982). It is known that the IBL growth is capped by the thermal inversion at the top of the ABL (Oke Reference Oke1976), 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 growth of the turbulence length scale $\hat {l}_f$ using the following relationship:

(5.17)\begin{equation} \hat{l}_f=\frac{\delta}{1+\delta/H}. \end{equation}

For small values of $x$, $\hat {l}_f\approx \delta$, while $\hat {l}_f\to H$ as the value of $x\to \infty$. Variation of $\hat {l}_f$ with $x$ based on (5.17) is shown in figure 10. For simplicity, we here assume that the value of $H$ is constant. It is, however, important to note that depending on the size of the wind farm and the level of thermal stratification in the inversion layer, the IBL may lead to the growth of the entire ABL, so in reality $H$ might change with $x$ (Allaerts & Meyers Reference Allaerts and Meyers2017; Wu & Porté-Agel Reference Wu and Porté-Agel2017).

Figure 10. Variation of velocity scale $\hat {u}_f$, length scale $\hat {l}_f$ and farm turbulent viscosity $\nu _{t,f}$ based on the modelling approached elaborated in § 5.5. Variation of $\nu _{t,f}$ for the LES data (case S0) is also shown. In the figure, $\hat {u}_f$ is normalised by $\sqrt {c_{ft}}$, $\nu _{t,f}$ is normalised by $\nu _{t,f}(x=x_1+L_f)$ and $\hat {l}_f$ is normalised by $H$. Vertical dotted lines show the location of turbine rows, and the vertical dashed line shows $x_1+L_f$ where the turbulent viscosity $\nu _{t,f}$ is maximum.

Next, we determine the turbulence velocity scale $\hat {u}_f$. Within the wind farm, turbulence is mainly generated due to the shear caused by turbine forcing. Therefore, inspired by (Calaf et al. Reference Calaf, Meneveau and Meyers2010), we use $\sqrt {c_{ft}}U_0$, where $U_0=1$, to estimate the turbulence velocity scale $\hat {u}_f$ 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 Reference Stieren and Stevens2022). The generated turbulence in the turbine wake usually peaks at around $5$ rotor diameters downstream (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009). The constant turbulence velocity scale $\hat {u}_f$ caused by turbine forcing is therefore assumed to be extended from $x_1$ to $x_N+5$ as shown in figure 10. Further downstream, $\hat {u}_f$ 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. (Reference Belcher, Jerram and Hunt2003), the velocity 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 $\hat {u}_f$ as

(5.18)\begin{equation} \hat{u}_f =\begin{cases} 0 & \text{if } x< x_1,\\ \sqrt{c_{ft}} & \text{if } x_1< x< x_1+L_f,\\ \sqrt{c_{ft}}\left(\dfrac{x-x_1}{L_f}\right)^{{-}1} & \text{if } x>x_1+L_f, \end{cases} \end{equation}

where $L_f=(x_N-x_1)+5$. Variations of $\hat {u}_f$ and $\nu _{t,f}$ are shown in figure 10, where the farm turbulent viscosity $\nu _{t,f}$ is given by

(5.19)\begin{equation} \nu_{t,f}=c_2\hat{u}_f\hat{l}_f, \end{equation}

with $c_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 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 $\nu _{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 for the wake of a single wind turbine (Scott et al. Reference Scott, Martínez-Tossas, Bossuyt, Hamilton and Cal2023).

5.6. Estimation of local deficit coefficient $\eta _n$

As discussed earlier, $\bar {u}_{h,n}=1-\eta _nU_d(x=x_n)$, where $\eta _n$ is the local deficit coefficient for WT$_n$s. A consequence of linearising the momentum equation (5.8) is that wake effects are linearly superposed (Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021). Therefore, one can write

(5.20)\begin{equation} \eta_{n} =\frac{\displaystyle\sum_{m=1}^{n-1} \bar{u}_{d,m}(x_n,y_n,z_h)}{\displaystyle\sum_{m=1}^{n-1} U_{d,m}(x_n)}, \end{equation}

where $\bar {u}_{d,m}(x,y,z)$ is the local velocity deficit at $(x,y,z)$ caused by WT$_m$s and $U_{d,m}(x)$ is already defined in (5.12). To find $\bar {u}_{d,mn}$, one can use a Gaussian distribution $C\exp (-r^2/2\sigma ^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 $\sigma$ is the characteristic wake width. Therefore, for a semi-infinite wind farm, we obtain

(5.21)\begin{equation} \bar{u}_{d,m}(x_n)=\lim_{K\to\infty}C_{mn}\sum_{k={-}K}^{K}\exp \left(-\frac{(y_n-y_m-ks)^2}{2\sigma_{mn}^2}\right), \end{equation}

where $C_{mn}$ is the maximum velocity deficit caused by each WT$_m$ at $x=x_n$, $k$ denotes the column number which varies from $-\infty$ to $\infty$ and $\sigma _{mn}$ is the width of WT$_m$ wakes at $x=x_n$. Solving (5.21) gives

(5.22)\begin{equation} \bar{u}_{d,m}(x_n)=\sqrt{2{\rm \pi}}C_{mn}\frac{\sigma_{mn}}{s_y}\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left({-}2{\rm \pi}^2\frac{\sigma_{mn}^2}{s_y^2}\right)}\right], \end{equation}

where $\vartheta _{3}[z,q]$ is the Jacobi theta function defined as $1+2\sum _{l=1}^{\infty }q^{l^2}\cos (2lz)$ (Whittaker & Watson Reference Whittaker and Watson2020). 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 $\vartheta _{3}[z\pm {\rm \pi},q]=\vartheta _{3}[z,q]$, so $y_n$ and $y_m$ in (5.22) can be the spanwise location of any turbines in WT$_n$ and WT$_m$, respectively.

By definition $U_{d,m}(x)=\langle \bar {u}\rangle _{d,m}(x)$, so we can write

(5.23)\begin{align} U_{d,m}(x_n)=\lim_{K\to\infty}\frac{C_{mn}}{Ms_y}\sum_{k={-}K}^{K}\int_{-\infty}^{+\infty}\exp \left(-\frac{(y-y_m-ks_y)^2}{2\sigma^2_{mn}}\right)\textrm{d}\kern0.05em y=\sqrt{2{\rm \pi}}C_{mn}\frac{\sigma_{mn}}{s_y}. \end{align}

Therefore, from (5.22) and (5.23), we conclude that

(5.24)\begin{equation} \bar{u}_{d,m}(x_n)=U_{d,m}(x_n)\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left({-}2{\rm \pi}^2\frac{\sigma_{mn}^2}{s_y^2}\right)}\right]. \end{equation}

The wake width $\sigma _{mn}$ can be simplified by $k^{*}(x_n-x_m)+\sigma _0$, where the wake expansion $k^{*}$ is typically around 0.02–0.04 for offshore conditions and $\sigma _0$ is around 0.2–0.3 (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014). For simplicity, we assume $k^{*}=0.025$ and $\sigma _0=0.25$ toreduce (5.24) to

(5.25)\begin{equation} \bar{u}_{d,m}(x_n)\approx U_{d,m}(x_n)\vartheta_{3}\left[{\frac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left(-\frac{(x_n-x_m+10)^2}{80s_y^2}\right)}\right]. \end{equation}

Therefore from (5.20) and (5.25), we conclude

(5.26)\begin{equation} \eta_{n} =\frac{\displaystyle\sum_{m=1}^{n-1} U_{d,m}(x_n)\vartheta_{3}\left[{\dfrac{{\rm \pi}(y_n-y_m)}{s_y}},{\exp \left(-\dfrac{(x_n-x_m+10)^2}{80s_y^2}\right)}\right]} {\displaystyle\sum_{m=1}^{n-1}U_{d,m}(x_n)}, \end{equation}

where $n>1$ (for $n=1$, $\eta$ is zero) and $U_{d,m}(x)$ is defined in (5.12). Figure 11 shows the variation of the local deficit coefficient $\eta _n$ with the row number $n$ for both the Aligned Baseline (A0) and the Staggered Baseline (B0) wind farms. It is interesting to note that the value of $\eta _n$ 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 $\eta$ 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 $\eta _n$ decays and approaches a constant value of 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 $\eta _n$ 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 $\eta _n$. The value of $\eta _n$ 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 11(b), we can see how different distributions of the local velocity coefficient $\eta _n$ lead to different distribution of local hub-height velocity $\bar {u}_{h,n}$. 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 11(b) is similar to the data reported in other works (e.g. Chamorro, Arndt & Sotiropoulos Reference Chamorro, Arndt and Sotiropoulos2011; Stieren & Stevens Reference Stieren and Stevens2022). 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 $\eta _n$ developed in (5.26), to any conceivable wind farm layout provided the layout fulfils the requirements defined in § 5.1.

Figure 11. Variation of (a) local deficit coefficient $\eta _n$ and (b) local hub-height velocity $\bar {u}_{h,n}$ with row number $n$ for the Staggered Baseline (S0) and the Aligned Baseline (A0) cases based on (5.26).

5.7. Estimation of shear term $C_x$ and veer term $C_y$

As discussed earlier in (5.10), $C_{x}=-\nu _{t,f} \partial ^2 U_0/\partial z^2$ and $C_{y}=-\nu _{t,f} \partial ^2 V_0/\partial z^2$. Based on (5.9), $C_{x}$ and $C_{y}$ can be written as

(5.27a)$$\begin{gather} C_{x}={-}\frac{\nu_{t,f}}{\nu_{t,0}}f_c(V_g-V_0), \end{gather}$$

(5.27b)$$\begin{gather}C_{y}=\frac{\nu_{t,f}}{\nu_{t,0}}f_c(U_g-U_0). \end{gather}$$

As a first-order approximation, the vertical changes in the streamwise and spanwise velocities across the ABL are proportional to $G\cos \theta _0$ and $G\sin \theta _0$, respectively. Here, $G=\sqrt {U_g^2+V_g^2}$ is the geostrophic wind speed, and the cross-isobar angle $\theta _0$ is the angle between the wind direction on the ground surface and the geostrophic wind direction. One may thus write

(5.28a)$$\begin{gather} U_g-U_0=c_3G\cos\theta_0, \end{gather}$$

(5.28b)$$\begin{gather}V_g-V_0=c_3 G\sin\theta_0, \end{gather}$$

where $c_3$ is a constant. Values of $G$ and $\theta _0$ can be obtained from the widely used geostrophic drag law, which relates surface properties (e.g. $z_0$ and $u_*$) to geostrophic wind speed on top of the ABL (e.g. Blackadar & Tennekes Reference Blackadar and Tennekes1968; Hess & Garratt Reference Hess and Garratt2002; van der Laan et al. Reference van der Laan, Kelly, Floors and Peña2020). For a neutrally stratified ABL, the geostrophic drag law reads as (Liu, Gadde & Stevens Reference Liu, Gadde and Stevens2021)

(5.29a)$$\begin{gather} A=\ln\left(\frac{u_*}{z_{0,0}|f_c|}\right)-\frac{\kappa G}{u_*}\cos\theta_0, \end{gather}$$

(5.29b)$$\begin{gather}B={\mp} \frac{\kappa G}{u_*}\sin\theta_0, \end{gather}$$

where $A$ and $B$ are universal empirical constants, and values of $A = 1.8$ and $B= 4.5$ used in the Wind Atlas Analysis and Application Program (Floors, Troen & Kelly Reference Floors, Troen and Kelly2018) are adopted here. The minus sign on the right-hand side of the second equality in (5.29) relates to the northern hemisphere where $f_c$ is positive, whereas the positive sign is for the southern hemisphere where $f_c$ is negative (Liu et al. Reference Liu, Gadde and Stevens2021). From (5.27), (5.28) and (5.29), we can therefore approximate the values of $C_x$ and $C_y$ as follows:

(5.30a)$$\begin{gather} C_{x}=c_3|f_c|\frac{\nu_{t,f}}{\nu_{t,0}}\frac{u_*}{\kappa}B, \end{gather}$$

(5.30b)$$\begin{gather}C_{y}=c_3f_c\frac{\nu_{t,f}}{\nu_{t,0}}\frac{u_*}{\kappa}\left[\ln\left(\frac{u_*}{z_{0,0}|f_c|}\right)-A\right]. \end{gather}$$

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