1. Introduction and background
Generating power from wind has gradually acquired an ever more important role in our society as a viable alternative source to fossil fuels. It is, therefore, increasingly important to gain a deeper understanding and more accurate predictability of the wind power availability, the aero-elastic fatigue loads on the wind turbine blades and the interaction of wake flow with atmospheric turbulence. It is well understood that a wind turbine, by extracting momentum from the wind, reduces the available power for the turbines downstream of it (Porté-Agel, Bastankhah & Shamsoddin Reference Porté-Agel, Bastankhah and Shamsoddin2020). Barthelmie et al. (Reference Barthelmie, Pryor, Frandsen, Hansen, Schepers, Rados, Schlez, Neubert, Jensen and Neckelmann2010) showed that this reduction can be of the order of 
$40\,\%$, presenting an industrial challenge when predicting the efficiency of a wind farm as a whole. Within the wake flow, the turbulent levels are enhanced due to the rotation of the blades and the presence of the turbine itself as a static obstacle to the incoming flow (Hancock & Pascheke Reference Hancock and Pascheke2014b; Graham Reference Graham2017; Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021). The problem is complicated by the fact that most of the work on wind power has considered winds in neutral conditions, where the (potential) temperature gradient across the atmospheric boundary layer (ABL) is absent; nevertheless, the ABL is rarely in neutral conditions. In the North Sea, for example, the winds are in neutral condition only approximately 30 % of the time (Barthelmie et al. Reference Barthelmie2011), and similar conditions are also typical of other offshore and onshore sites (Floors et al. Reference Floors, Gryning, Peña and Batchvarova2011; Wharton & Lundquist Reference Wharton and Lundquist2012; Peña et al. Reference Peña, Floors, Sathe, Gryning, Wagner, Courtney, Larsén, Hahmann and Hasager2016; Ichenial, El-Hajjaji & Khamlichi Reference Ichenial, El-Hajjaji and Khamlichi2021). The effect of thermal stability has been shown in physical (laboratory and full-scale) and numerical experiments to affect the wake recovery rate within the wake flow of single and multiple turbines (Abkar & Porté-Agel Reference Abkar and Porté-Agel2013; Tabib, Rasheed & Kvamsdal Reference Tabib, Rasheed and Kvamsdal2015; Wu & Porté-Agel Reference Wu and Porté-Agel2017; Allaerts & Meyers Reference Allaerts and Meyers2018; Abkar, Sørensen & Porté-Agel Reference Abkar, Sørensen and Porté-Agel2018). The consensus resides on the fact that the velocity deficit in the wake flow assumes an exponential decay modulated by the thermal stability, although a less convincing agreement is reached on the parameters describing such decay (e.g. Magnusson & Smedman Reference Magnusson and Smedman1994; Barthelmie et al. Reference Barthelmie2004; Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2010; Zhang, Markfort & Porté-Agel Reference Zhang, Markfort and Porté-Agel2012; Hancock & Pascheke Reference Hancock and Pascheke2014b; Abkar & Porté-Agel Reference Abkar and Porté-Agel2015; Machefaux et al. Reference Machefaux, Larsen, Koblitz, Troldborg, Kelly, Chougule, Hansen and Rodrigo2016; Stevens & Meneveau Reference Stevens and Meneveau2016; Bodini, Zardi & Lundquist Reference Bodini, Zardi and Lundquist2017; Xie & Archer Reference Xie and Archer2017; Ning & Wan Reference Ning and Wan2019; Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021; Gadde & Stevens Reference Gadde and Stevens2021).
Less explored in the literature are higher-order turbulence quantities and how these are affected by thermal stability. Cal et al. (Reference Cal, Lebrón, Castillo, Kang and Meneveau2010) investigated a relatively large (
$3\times 3$) array of turbines, considering turbulence properties in the near field at several locations within the array. They concluded that the vertical fluxes of mean kinetic energy are important if one wants to estimate the overall power extraction. Further analysis of the same dataset considered the most relevant fluxes of kinetic energy across a wind turbine stream tube, confirming the key role of turbulence fluxes in shaping its dynamics (Lebron, Castillo & Meneveau Reference Lebron, Castillo and Meneveau2012). Recent attempts applied a novel methodology to identify characteristic regions of the flow within a large wind farm based on the structure of the Reynolds stress anisotropy. Through a clustering technique, Ali et al. (Reference Ali, Hamilton, Calaf and Cal2019a) argued that major perturbations control the operation of the wind farm during the diurnal cycle, with the effect of thermal stability governing the exchange of momentum flux. Chamorro & Porté-Agel (Reference Chamorro and Porté-Agel2009) focused on characterising the turbulence intensity in the wake of a single wind turbine, confirming that the presence of the latter vastly enhances the turbulence levels in the upper part of the wake. In the lower part of the wake instead, turbulence intensities are reduced when compared to the undisturbed flow. By investigating spectra at different wall-normal heights in the wake flow, they also correlated the turbulent intensity modification due to the turbine with relatively large scales in the flow. They extended this work to a large laboratory-scale wind farm, investigating similar quantities and overall power production (Chamorro, Arndt & Sotiropoulos Reference Chamorro, Arndt and Sotiropoulos2011). Chamorro et al. (Reference Chamorro, Guala, Arndt and Sotiropoulos2012) introduced the concept of a wind turbine as an ‘active filter’, where it was found to energise the high-frequencies, while also disrupting the large and very-large scales of the flow in proximity of the wall. Fuertes, Markfort & Porté-Agel (Reference Fuertes, Markfort and Porté-Agel2018) reported results from full-scale wind turbine wakes invested by different wind conditions. They found that higher incoming turbulence levels help the entrainment of external flow and mixing resulting in faster wake recovery. Gambuzza & Ganapathisubramani (Reference Gambuzza and Ganapathisubramani2023) decoupled the effect of the free stream turbulence from the integral time scales on the wake flows, showing that higher integral time scales result in a slower wake evolution (mimicking lower-turbulence inflow conditions) irrespective of the turbulence ambient levels.
The topic of turbulent correlations in wake flows, however, is vastly unexplored in the current literature, particularly as a function of the atmospheric state, and it is absent from recent reviews (Stevens & Meneveau Reference Stevens and Meneveau2016; Porté-Agel et al. Reference Porté-Agel, Bastankhah and Shamsoddin2020). These correlations, however, have obvious repercussions on the unsteady turbulence loads experienced by the turbines (Dörenkämper et al. Reference Dörenkämper, Witha, Steinfeld, Heinemann and Kühn2015; Tian, Ozbay & Hu Reference Tian, Ozbay and Hu2018) and their life cycle, and are, therefore, of great industrial interest. Meng et al. (Reference Meng, Lien, Glinka, Li and Zhang2018) pointed out that wake-induced fatigue can be responsible for up to a 15 % increase in the rotor loads. Liu & Stevens (Reference Liu and Stevens2021) examined the combined effects of thermal stability and complex terrain on power production, confirming that different stability conditions result in varied fatigue loads. They also highlight how thermal characteristics and complex terrain effects can coexist, giving rise to non-trivial interactions. Understanding velocity correlation maps is also important in practical wind resource assessment methodologies as full wind farm solvers still require prohibiting computing resources (Stevens & Meneveau Reference Stevens and Meneveau2016; Cheng & Porté-Agel Reference Cheng and Porté-Agel2018; Ichenial et al. Reference Ichenial, El-Hajjaji and Khamlichi2021). As an example, the methodology proposed by Veers (Reference Veers1984) is a simple – yet elegant – solution that relies on simulating the instantaneous three-dimensional wind field impacting on a wind turbine based solely on information from the power spectral density (PSD) of the incoming wind and the two-point velocity correlations in space across the turbine diameter (i.e. the ‘wind coherence function’). It was shown that considering the rather simplistic assumptions of the model and its relatively low computational cost, which renders it very appealing for industrial applications, numerical predictions agreed surprisingly well with the field measurements over a large sample of test cases (Veers Reference Veers1984, Reference Veers1988). Experimental work on the topic of velocity correlations in stratified flows is sparse, but not completely unexplored. Howard, Chamorro & Guala (Reference Howard, Chamorro and Guala2016), most notably, considered different thermal regimes (from weakly stable to weakly convective), focusing predominantly on the spatial evolution of the large-scale motions confirming that, in the very near field ahead of a wind turbine, the inclination of the streamwise velocity correlation maps is linked to the thermal characteristic of the incoming boundary layer. Their investigation, however, is limited to the near field and does not consider either the effect of increasingly stratified flows or the lateral coherency of the turbulence. Ali et al. (Reference Ali, Hamilton, Calaf and Cal2019b) numerically investigated the effect of thermal stratification on the turbulent kinetic energy budget, and examined triple velocity correlations (i.e. the turbulent convection) and their role in energising the wake flow. For stable cases, the variation in the correlation was found to be vastly negligible for the wall-normal and spanwise components, while for convective winds, all three components of the velocity were deemed important to capture the energy budget. Also noteworthy is the work of Ali et al. (Reference Ali, Cortina, Hamilton, Calaf and Cal2017), which applying proper orthogonal decomposition, highlighted how the shape and size of the key energetic modes (a proxy for the turbulence length scales) are affected by the thermal state of the wind. Unstable and neutral cases were found to be dominated by large-scale coherent motions linked to the boundary layer scales, which were beneficial for the wake recovery. In contrast, for stably stratified flows, modes linked to the presence of the wind turbine acquired more importance.
Having considered the shortcomings of the existing literature, the aim of the current work is twofold. First (and foremost), it seeks to undertake a systematic investigation of three-dimensional velocity correlation maps in the near- and far-wake flows of wind turbines in industrially relevant stably stratified boundary layers. Second, given the lack of any theoretical restrictions for the application of Veers’ methodology (or similar) to these flows, it aims to highlight some necessary considerations for the future application of these methodologies to non-neutral winds. Laboratory data acquired in the presence of single and multiple wind turbines (up to two) and with increasing levels of thermal stability are presented herein to shed light on their combined effects on the wake flow characteristics. The investigation by Hancock & Pascheke (Reference Hancock and Pascheke2014b), carried out at the same facility, only examined single-point measurements for a single wind turbine. In addition, several improvements in the modelling technique have been made since then, including those by Hancock & Hayden (Reference Hancock and Hayden2018, Reference Hancock and Hayden2020, Reference Hancock and Hayden2021).
2. Experimental facility, instrumentation and details
2.1. Facility and boundary layer development
Experiments were conducted in the EnFlo wind tunnel at the University of Surrey. The facility is a suck-down tunnel with a bell-mouth inlet, open return and a large working section measuring 
$20\ {\rm m}\times 3.5 \ {\rm m} \times 1.5\ {\rm m}$ (
${\rm length}\times {\rm width}\times {\rm height}$). The uniqueness of this facility is that it enables simulations of atmospheric boundary layers under neutral, stable and convective conditions thanks to the inclusion of heating and cooling elements on its floor (and sidewalls), as well as the possibility of specifying a prescribed temperature profile on the incoming wind via inlet heaters. Further details on the facility are provided by Hancock & Hayden (Reference Hancock and Hayden2018) together with a comparison of the typical boundary layer development with customary scaling laws adopted in meteorology. Artificially thickened boundary layers are obtained by the use of 13 truncated triangular flat plate spires located at the tunnel inlet (
$x=-10$ m). These triangles have dimensions of 60 mm at the base, 4 mm at their apex and stand 600 mm tall, whilst being spaced 266 mm in the spanwise direction. To recreate realistic surface roughness conditions typical of offshore seas (Counihan Reference Counihan1975), roughness elements 50 mm wide, 16 mm high and 5 mm thick were employed in a staggered arrangement and covered the entire tunnel fetch, starting 10 m upstream of the first turbine. This boundary layer set-up was identical to that described by Hancock & Hayden (Reference Hancock and Hayden2018, Reference Hancock and Hayden2021) (with a high degree of repeatability in the mean and high-order quantities) and resulted in the roughness lengths reported in table 1.
Table 1. Summary of boundary layer characteristics at 
$x/D=0$. Free stream inversion gradient is 
$20\ {\rm km}^{-1}$ for both stable cases. Here, 
$u_*$ is the friction velocity, 
$\theta ^*$ is the temperature scale in 
$K$, 
$z_0$, 
$z_{0\theta }$ and 
$L_0$ are the surface-layer roughnesses and surface Obukhov length in mm, 
$Re_\delta$ is the boundary layer thickness Reynolds number, and 
$(\overline {w\theta })_0$ is the kinematic heat flux at the surface in 
$\textrm{Km\ s}^{-1}$. Thermal quantities are taken from Hancock & Hayden (Reference Hancock and Hayden2018, Reference Hancock and Hayden2021).
2.2. Stable boundary layers: scaling consideration and implementation
The conventional scaling between full- and laboratory-scale boundary layers is to match the ratios between the Obukhov length and the turbine rotor diameter 
$L/D$ (assuming that the hub height is also a function of the rotor diameter), and that of the boundary layer depth to the surface Obukhov length 
$\delta /L_0$. The first is defined as
\begin{equation} \frac{L}{D}={-}\frac{1}{\kappa} \frac{\bar{T}}{g} \frac{(-\overline{u w})^{3 / 2}}{\overline{w \theta}} \frac{1}{D}, \end{equation}
where 
$\bar {T}$ is the mean absolute temperature, 
$\theta$ is the fluctuating potential temperature, 
$\overline {w \theta }$ is the kinematic vertical heat flux, 
$g$ is the acceleration due to gravity, and 
$\kappa$ and 
$-\overline {u w}$ are the von Kármán constant and the Reynolds shear stress, respectively. At the wall, 
$L=L_{0}$, so that
\begin{equation} \frac{L_{0}}{D}={-}\frac{1}{\kappa} \frac{\bar{T}}{g} \frac{u_{*}^{3}}{(\overline{w \theta})_{0}} \frac{1}{D}, \end{equation}
where 
$u_{*}$ is the friction velocity and 
$(\overline {w \theta })_{0}$ denotes the surface kinematic heat flux. Assuming that the temperature fluctuations are driven by the mean temperature gradient, Hancock & Pascheke (Reference Hancock and Pascheke2014a) have shown that the latter (in the mixed layer and overlying inversion) for full-to-model scale similarity must satisfy
\begin{equation} \left(\frac{D}{U_R}\right)^2\frac{\partial\varTheta}{\partial z}=\mathrm{const.}, \end{equation}
where 
$U_R$ is a reference velocity and 
$\varTheta$ is the mean potential temperature. Hence, as 
$D/U_R$ is much smaller in the laboratory than at full scale, it follows that the temperature gradient must be increased according to (2.3).
This work considers three thermal stability conditions: neutral, stable without an overlying inversion and stable with an overlying inversion. To generate stable boundary layers, the scaling above is applied coupled with the methodology of Hancock & Hayden (Reference Hancock and Hayden2018, Reference Hancock and Hayden2020). This implies matching the spires and the wall roughness of the neutral case, whilst also heating the flow at the working-section inlet to a prescribed temperature profile and cooling a fetch of the floor to a fixed temperature. The cooled fetch, here, spanned in the range 
$6< x<19$ m. It should be noted that an initial fetch of the tunnel floor is kept uncooled; this is necessary to obtain smoothly varying turbulence statistics across the boundary layer depth (Hancock & Pascheke Reference Hancock and Pascheke2014a, Reference Hancock and Paschekeb). It can be shown how this arrangement generates well-behaved atmospheric boundary layers respecting the Monin–Obukhov surface-layer similarity (Hancock & Pascheke Reference Hancock and Pascheke2014a; Hancock & Hayden Reference Hancock and Hayden2018, Reference Hancock and Hayden2020). The inlet temperature profile for the inversion case is prescribed based on the stable case profile and superimposing an additional linear increase of 2 K per each wall-normal interval at the inlet heaters (
$\Delta z=100$ mm). This superposition started at height 
$z=250$ mm so that the inversion gradient above the boundary layer is 
$20\ {\rm K}\ {\rm m}^{-1}$ (see discussion in § 3.1). The thermal boundary layer characteristics are reported in table 1.
2.3. Instrumentation
Single-component velocity measurements and two-point velocity correlations were acquired with two recently calibrated laser Doppler anemometers (FiberFlow, Dantec Dynamics, Denmark) at a minimum frequency of 100 Hz for 120 s for each measurement point. Tests were also conducted with longer acquisition time, however, these showed that 120 s were sufficient to ensure statistical convergence of the data. Statistical errors are typically within 
$\pm 0.5\,\%$ and 
$\pm 5\,\%$ for the mean and second-order moments within a 95 % confidence level. The laser Doppler anemometer (LDA) only measured the velocity component 
$U$ along the streamwise direction 
$x$. Spanwise and wall-normal components (along 
$y$ and 
$z$, respectively) were deemed to be of secondary importance herein as per the hypotheses of Veers (Reference Veers1988). Mean velocities and turbulent fluctuations are indicated by capital and lower-case symbols, respectively, while time-averaging operations are indicated by an overbar. For convenience, we employ a relative wall-normal position, 
$Z$, with respect to the turbine's hub, where 
$Z=z-z_h$, so that 
$Z=0$ is the hub height. The two LDA probes (with a focal length of 
$300$ mm) were mounted on two independently controlled traverse systems with sub-mm accuracy operating on carriages running internally along the roof of the wind tunnel. The wind-tunnel free stream velocity, 
$U_\infty$, was set to 
$1.5\ {\rm ms}^{-1}$, as measured via means of an ultrasonic anemometer mounted at 
$(x,y,z)=(-5\ {\rm m},0.7\ {\rm m},1\ {\rm m})$. The origin of the right-hand coordinate system is positioned 10 m downstream from the inlet of the working section, at hub height on the wind tunnel centreline. Due to the tunnel design (i.e. absence of settling chamber, contraction and screens), the flow uniformity across a spanwise distance of 1 m is within 
$0.7\,\%$ and 
$7\,\%$ for 
$U$ and 
$u^2$, respectively. Measurements are automated and coordinated via in-house software (LabView, National Instruments, United Kingdom).
2.4. Two-point velocity correlations and turbulence length scales
Two LDA probes were employed to measure the two-point correlation of streamwise fluctuating velocity components. Here, we use the customary definition of the normalised correlation coefficient as
\begin{equation} \rho=\frac{\overline{u_1u_2}}{(\overline{{u_1}^2}\ \overline{{u_2}^2})^{1/2}}. \end{equation}
Once the correlation coefficient is evaluated, this can be used as a proxy to define a turbulence length scale, defined as the spacing between points characterised by 
$\rho = 0.5$, similarly to Christensen & Wu (Reference Christensen and Wu2005). We can therefore define streamwise, spanwise and wall-normal length scales, namely 
$L_{x}$, 
$L_{y}$ and 
$L_{z}$, based on the correspondent correlation coefficients.
2.5. Wind turbine models
The wind turbine models used in this work were developed by Hancock & Pascheke (Reference Hancock and Pascheke2014b). The turbine is a 1 : 300 scale representation of a 5 MW wind turbine typical of off-shore applications. The radius of the turbine is 
$R=208$ mm, with the hub at 
$z_h=300$ mm, where the 
$_h$ subscript indicates hub quantities. Therefore, the turbine rotor spans 
$92\ {\rm mm} < z<508\ {\rm mm}$. The blades are tapered and twisted 0.8 mm thick flat plates manufactured in-house out of carbon fibre with layout (
$45^\circ, 0^\circ, 0^\circ, 45^\circ$) to the blade radial axis. The turbine and blades geometries are shown in figure 1. The rationale for this flat plate-like geometry is to account for the blades’ aerodynamics at the operating laboratory Reynolds number, much lower than that of the full-size machine (Hancock & Pascheke Reference Hancock and Pascheke2014b). The pitch of the blades is adjustable and it is specified in figure 1(b). This was measured across all blades and all turbines to be within 
$\pm 0.5^\circ$ of the design angle at the beginning of the tests, with the system described by Farr (Reference Farr2014). An example of the measured data is offered in figure 1(b). The design tip-speed ratio 
$TSR=\varOmega R/U_h=6\pm 1.5\,\%$, throughout the tests, where 
$\varOmega$ is the angular velocity and 
$U_h$ is the mean velocity at hub height. For the downstream turbine, this velocity was obtained by area-averaging the time-averaged velocity in the wake flow across the rotor diameter measured by the LDA system at the beginning of the experiment. To achieve this 
$TSR$, a 4 : 1 gear was mounted onto the motor's shaft. With these settings, an estimated thrust coefficient of 0.48 was achieved. The blades were mounted on a 13 mm diameter hub, itself mounted directly onto the motor and gearbox shaft with the same diameter and length of 60 mm. A motor (RE-13-118617, 6V, 3W, Maxon, Switzerland) in combination with a four-quadrant controller (4-Q-DC-Servoamplifier LSC, 0-30W, Maxon, Switzerland) and digital encoder was used to control the turbines’ rotation. The sense of rotation was clockwise when viewing the turbines from upstream. The drivetrain system was secured with an unex clip to the tower of circular cross-section, also 13 mm in diameter. The tower was connected to a small rectangular base secured to the wind tunnel floor.
Figure 1. (a) Wind turbine sketch with overall dimension (not to scale). (b) Chord and blade twist characteristics across a single blade.
2.6. Test cases and further measurements details
Three different turbine configurations were considered in this work, alongside three thermal stability conditions, for a total of nine cases, as in table 1. The case named 0WT represents the baseline case (i.e. the undisturbed boundary layer). Cases 1WT and 2WT considered one and two wind turbines, respectively. These are positioned in line at 
$x=10\,000$ mm and 
$x=12\,496$ mm (i.e. 
$x=0D$ and 
$x=6D$). A sketch of the experimental set-up is provided in figure 2. Two-point synchronised velocity measurements were acquired with the reference probe situated at different wall-normal and spanwise locations (
$z=z_h-0.35D, z_h, z_h+0.35D$ and 
$y=-0.35D, 0, 0.35D$, respectively). For each case, streamwise locations from 
$x=0D$ to 
$x=15D$ were considered. The exact measurement stations are reported in table 2, with the same symbols but different colours (red and blue) indicating comparable stations (e.g. 
$x=1D$ for the 1WT case and 
$x=7D$ for the 2WT are acquired at the same relative distance from respective turbines).
Figure 2. Schematic of the experimental set-up (not to scale). The red dot indicates the origin of the coordinate system, and the reference anemometer sets the free stream speed.
Table 2. Summary of data and symbols. 0WT will be in black, 1WT will be in red and 2WT will be in blue. Solid symbols indicate the location upstream of the relevant turbine, while empty symbols indicate wake-flow measurements. Colour shades in figure 4(c) indicate measurements not at hub height. Neutral, stable without an overlying inversion and stable with an overlying inversion stability conditions are presented as solid (—), dashed (
$--$) and dotted (
$\cdots$) lines and labelled (N), (S) and (MI), respectively.
3. Results and discussion
3.1. Effect of thermal stability on the undisturbed boundary layer
This section focusses on the changes to the undisturbed boundary layer due to different thermal conditions. First, we discuss mean and fluctuating characteristics in § 3.1.1, then the analysis considers velocity correlations in § 3.1.2.
3.1.1. Mean profiles and turbulence statistics
Baseline boundary layer measurements at different streamwise stations in the absence of wind turbines and thermal stability are presented in figure 3(a). The collapse of all curves across the streamwise stations demonstrates that the boundary layer is fully developed (free stream velocity change across the measurements domain is below 
$2\pm 0.7\,\%$). Furthermore, the mean velocity is sufficiently homogeneous (variation across the span below 5 %) with profiles of Reynolds stresses and heat flux also horizontally homogeneous over a large streamwise fetch (omitted here but reported by Hancock & Hayden Reference Hancock and Hayden2018, Reference Hancock and Hayden2021). The effect of thermal stability (mean temperature profiles reported in figure 3d) on the mean velocity profiles are presented in figure 3(b) and normalised by the outermost local velocity measured, 
$U_e$. This is employed (due to limitations in the measurement set-up) in lieu of the local free stream velocity by Hancock & Hayden (Reference Hancock and Hayden2021); these are, however, within 2 %. There is no more than a slight change in the mean velocity profiles with the introduction of stable conditions, even in the presence of an overlying inversion. This is partly due to the use of spires that somewhat prescribe the depth of the boundary layer irrespective of its thermal characteristics. The effect of the different thermal states becomes obvious when considering the turbulent fluctuations in figure 3(c). Here, it is clear how increasing the thermal stability suppresses the turbulent fluctuations across the entire wall-normal range, in line with previous literature (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017). Of interest is the concurrence of the fluctuating profiles for the two stable cases (dashed and dotted lines) below 
$Z/D<-0.4$. This has been observed previously by Hancock & Hayden (Reference Hancock and Hayden2020), who remarked on how changing the imposed inversion condition left the bottom third of the boundary layer almost unaffected.
Figure 3. (a) Mean streamwise velocity development in neutral conditions. (b) Mean and (c) fluctuations in different thermal conditions. (d) Imposed temperature profiles. Power spectral density of streamwise velocity fluctuations for (e) 
$Z=-0.5D$ and ( f) 
$Z=0.5D$. All data are taken at 
$y/D=0$. The spectral slope represents the inertial subrange (
$\beta =-5/3$).
The effect of thermal stability on the power spectral density of the velocity fluctuations for two different heights at 
$x/D=1$ are shown in figure 3(e, f). These heights are chosen to aid the discussion in § 3.2.4. The two stable cases in figure 3(e) possess lower energy across the entire resolved frequency range; this is due to the attenuation of turbulent fluctuations associated with thermal stability just discussed (figure 3c). Of interest is the very mild effect of the temperature inversion at this wall-normal height (compare dashed and dotted lines). A similar effect is shown at hub height (omitted here), but in figure 3( f), the effect of the temperature inversion becomes more significant. Overall, the larger effects are seen in the lowest frequencies and highest boundary layer position; an indication that the thermal stratification involves the large and very large scales in the flow (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017; Hancock & Hayden Reference Hancock and Hayden2020). This effect of stability in dampening turbulence fluctuations and its bias towards large scales (or low frequency) have repercussions on the overall power that can be extracted from the wind; first, as the power a turbine can produce is linked to the turbulence intensity characterising the flow (Elliott & Cadogan Reference Elliott and Cadogan1990) and, second, given that wind turbines are known to be more effective at harvesting power from velocity fluctuations at low frequencies (Sheinman & Rosen Reference Sheinman and Rosen1992; Chamorro et al. Reference Chamorro, Lee, Olsen, Milliren, Marr, Arndt and Sotiropoulos2015; Gambuzza & Ganapathisubramani Reference Gambuzza and Ganapathisubramani2021).
3.1.2. Velocity correlations
Before the effect of thermal stability on two-point velocity correlations is discussed, it is useful to introduce these quantities in neutral boundary layers. All vertical correlations are presented at three reference heights: 
$z_{ref}=z_h-0.35D$ (
$Z_{ref}=-0.35D$), 
$z=z_h$ (
$Z_{ref}=0$) and 
$z=z_h+0.35D$ (
$Z_{ref}=0.35D$). The streamwise, and lateral directions are also normalised with the turbine diameter to facilitate comparisons.
The vertical correlations in figure 4(a) at the different wall-normal reference heights present a correlation coefficient 
$\rho$ that decays (although remaining positive) nearly symmetrically with the wall-normal distance from its value of unity at zero separation (i.e. 
$z=z_{ref}$). The only notable difference between the heights is that the coherence profile gets broader as the reference height increases; this is consistent with rough wall turbulent boundary layers theory that prescribes an increase of the characteristic turbulence length scales with wall-normal distance from the wall, following the attached-eddy hypothesis (Sillero, Jimenez & Moser Reference Sillero, Jimenez and Moser2014; Placidi Reference Placidi2015; Marusic & Monty Reference Marusic and Monty2019). To quantify it, a wall-normal correlation length scale, 
$L_{z}$, can be defined as in § 2.4. In figure 4(a), this is found to be 
$L_{z}/D=0.30,0.35,0.36$ for 
$z=z_h-0.35D, z_h, z_h+0.35D$. This length scale is found to be largely independent of 
$x$ once the boundary layer is fully developed and in equilibrium with the underlying surface (see concurrence of data for 
$x=1D, x=6D$ and 
$x=15D$ in figure 4a); this is also in accordance with the literature (Adrian, Balachandar & Lin Reference Adrian, Balachandar and Lin2001; Christensen & Adrian Reference Christensen and Adrian2001; Volino, Schultz & Flack Reference Volino, Schultz and Flack2007). It must be pointed out that one of the features of this work is the fact that the boundary layer thickness does not vary considerably with the thermal characteristics of the wind. This can have repercussions on the vertical length scales of the turbulence.
Figure 4. (a) Vertical correlation coefficient for the three reference heights: 
$Z_{ref}=-0.35D, 0, 0.35D$. (b) Lateral and (c) longitudinal correlation coefficients at 
$y_{ref}=0$ and 
$Z_{ref}=0$. Data points in panel (b) are mirrored around 
$y/D=0$. (c) 
$Z_{ref}=-0.35D$ (dark grey) and 
$Z_{ref}=0.35D$ (light grey). Dashed lines indicate 
$\rho =0$.
The lateral correlation coefficients in figure 4(b) present a similarly sharp decrease from their maximum, reaching negative values for 
$0.4\lessapprox |y/D|\lessapprox 2$; this negative correlation is a consequence of continuity. Regardless of the streamwise location considered, a concurrence of the statistics can be observed, with 
$L_y/D\approx 0.32$. This is comparable to the average vertical length scale in agreement with canonical boundary layers. Finally, the longitudinal correlation coefficient is shown in figure 4(c). Its associated averaged length scale across the streamwise stations is the largest at 
$L_x/D=0.94$ for 
$z=z_h$. This correlation is related to the size of the large scales permeating the boundary layer (Adrian, Meinhart & Tomkins Reference Adrian, Meinhart and Tomkins2000; Adrian & Marusic Reference Adrian and Marusic2012), hence it is of no surprise that this is found to be approximately three times larger than 
$L_z$ and 
$L_y$ at 
$z=z_h$. It is also noticeable how 
$L_x$ is affected by the reference height. This is a well-known feature of rough-wall turbulence (Volino, Schultz & Flack Reference Volino, Schultz and Flack2009; Placidi Reference Placidi2015).
Now that the general features of turbulence length scales have been discussed for the baseline neutral case, figure 5 shows the effect of the thermal stability on these length scales. The stability affects the vertical correlation in figure 5(a), where an increase in stability is accompanied by a more rapid reduction in the correlation coefficient with height. It is important to note how this coefficient in the bottom part of the boundary layer is fairly unaffected by modification of the inversion conditions (see concurrence of dashed and dotted lines centred around 
$Z/D=-0.4$). The findings discussed so far extend the observations by Hancock & Hayden (Reference Hancock and Hayden2020), who found that a change of imposed inversion condition left the (mean and turbulence) profiles in the bottom third of the boundary layer largely unaffected to velocity correlations and spectra. The top of the boundary layer (
$Z/D>0$), however, is significantly influenced by the inversion; the reduction in 
$L_z$ is 26 %, and 51 % compared to the neutral case, at hub height and 
$Z=0.35D$, respectively. Even without the inversion, this length scale is reduced by 11 % and 17 % by the thermal effects at the same heights. As for the findings in figure 4, these correlation coefficients were found to be independent of 
$x$ (omitted here). Figure 5(b) shows the sensitivity of the lateral correlation coefficient to the thermal stability; this is near null for 
$-0.5D\leqslant y\le 0.5D$, and changes are small outside this range. Figure 5(c) presents the longitudinal correlation coefficients. Once again, a higher stability results in a lower longitudinal coherency. Reductions in the longitudinal length scale, 
$L_x$, were found to be 23 % and 43 % without and with the overlying inversion, compared to the non-stratified case. Similar reductions in the correlation are observed at all wall-normal heights, although these are omitted for clarity from figure 5(c).
Figure 5. (a) Vertical, (b) lateral and (c) longitudinal correlation coefficients. Data points in panel (a) are for the three reference heights: 
$Z_{ref}=-0.35D,0,0.35D$, whilst those in panel (b) are mirrored around 
$y/D=0$. Measurements in panels (b,c) are taken at hub height. Dashed lines indicate 
$\rho =0$. (d) A subset of the data in panel (a) with the inclusion of error bars.
Finally, a note on the measurement uncertainty and noise. Throughout the manuscript, we will only discuss changes across cases when results show discernible trends as a function of the thermal stability or due to the presence of the wind turbines despite uncertainty and noise. An example of the typical statistical uncertainty is provided via error bars in figure 5(d), which is a subset of the data presented in figure 5(a). It is clear how the previously discussed trends still hold irrespective of the measurement uncertainty and the signal-to-noise ratio. Hereafter, error bars are omitted, however, we employ markers whose size is conservatively indicative of the uncertainty in the measured quantities.
3.2. Wake flows
Having discussed the overall effect of thermal stability on the undisturbed boundary layer in § 3.1, we now shift our attention to the wake flows, where measurements for a single turbine are shown in red and those for two turbines in blue.
3.2.1. Velocity correlations in neutral flow
We proceed with the analysis of two-point correlation coefficients (and associated length scales) in the wake. Figure 6 shows the effect of wind turbines on the vertical profiles of the streamwise velocity correlation coefficients at the same reference heights of figure 4(a). To facilitate the discussion, these are presented separately and undisturbed flow data, at 
$x/D = 6$, are also provided (however, this particular location is irrelevant given the concurrence of the statistics shown in figure 3). The first point to note is that the coefficient upstream of the single turbine is largely unaffected by its presence, as can be seen from the concurrence of the * symbols with the black lines (at all heights). This is perhaps surprising given that the mean velocity profiles upstream, but close to the turbine, are affected by the impedance of the turbine (see figure 14(a) and discussions of Hancock & Farr Reference Hancock and Farr2014; Graham Reference Graham2017). A clear feature of the coefficients in figure 6 is their significant narrowing downstream of the turbines (see locations 
$x/D=2,5.8$). This is due to the disruption of the undisturbed turbulence given by the wake of the rotor and it is consistent with the smaller integral length scale of the wake-flow turbulence (España et al. Reference España, Aubrun, Loyer and Devinant2012; Hancock & Farr Reference Hancock and Farr2014; Tang et al. Reference Tang, Lam, Shum and Li2019). For example, the 
$L_{z}$ length scale at hub height (in figure 6b) is reduced by 50 % 
$5.8D$ downstream of the first turbine. It is also noticeable from figure 6(a) how in the near wake, 
$x/D=2$, the correlation becomes slightly negative due to the characteristic annular shape of the wake; this feature is lost farther downstream, 
$x/D=5.8$, as the wake becomes more homogeneous across the turbine rotor, but is also seen some six diameters downstream of the second turbine.
Figure 6. Vertical correlation coefficients for (a) 
$Z_{ref}=-0.35D$, (b) 
$Z_{ref}=0$ and (c) 
$Z_{ref}=0.35D$. Grey lines represent the dimensions of the rotor, while the dashed lines indicate 
$\rho =0$.
Figure 6(b) shows that the effect of the second turbine on the upstream flow is minimal from the closeness of the two sets of measurements for 
$x/D = 5.8$. However, there is a curious contrast in the profiles at a matched relative distance from their closest turbine (
${\bigtriangleup}$ symbols in red and blue). Below hub height, these are nearly concurrent, however, above it, the correlation coefficient in the case of 2WT has returned to that of the undisturbed flow, despite changes still present in figures 14(c) and 14(d); this is unexpected and perhaps fortuitous. Even more marked is the behaviour in figure 6(c), where the profile at 
$x/D = 11.8$ closely resembles the undisturbed profile in black, both above and below the reference height.
It can be concluded that the vertical correlation is reduced by the presence of the first wind turbine, but its recovery can potentially be enhanced by the added turbulence from the second turbine. Alternatively, this could indicate that the effect of the second turbine is not as strong as the first one due to the wake flow it is operating within. It would be interesting to investigate the correlation shape further downstream of the wake-impinged turbine. More widely, there is also an interesting question over what would be the shape of the correlation coefficient curves behind third and subsequent turbines. Overall, the recovery to undisturbed flow conditions is faster in the outer part of the boundary layer due to the elevated levels of turbulence.
Let us now consider the lateral correlation coefficients in figure 7(a). In the wake, these profiles are not symmetrical because of the wake rotation. The presence of the first turbine reduces the lateral coherency of the turbulence both upstream and downstream. The 
$L_{y}$ length scale at hub height is reduced by nearly 40 % downstream of the first turbine (
$x=5.8D$), compared to the reduction of 50 % in 
$L_z$ at the same location. The correlation measured around the second turbine is broader than that of the undisturbed case for 
$y/D>0$ and narrower for 
$y/D<0$ (with 
$L_y$ overall only reduced by 2 % compared to the undisturbed case); this is particularly evident in the correlation tails for 
$y/D>0$. Given the sense of rotation discussed in § 2.5, 
$y<0$ corresponds to a rotor-induced positive wall-normal velocity component so that the blades lift up low-momentum fluid away from the wall. The asymmetry is much stronger after the second turbine which is consistent with a superposition of effects from the turbines rotating in the same direction. The asymmetry is absent upstream of the first turbine (*-profile in red for 
$x/D = -0.2$), which is consistent with blade-element theory, prescribing no flow-induced rotation upstream of the rotor. The lateral broadening of the correlation around the second turbine discussed above is particularly interesting as to why it arises. Here, the wake features dominate the problem, allowing the flow not only to match the correlation coefficient typical of the undisturbed flow but also to exceed it; this is rather curious and it is not clear what would happen further downstream. When it comes to the upstream effect of the second turbine (blue *-profile for 
$x/D = 5.8$), there is again a clear upstream influence and it is in contrast to that seen in the vertical profiles, where no effect was seen. This time, the influence is asymmetric, with a pronounced negative lobe on the 
$y$-negative side, but not on the 
$y$-positive side.
Figure 7. (a) Lateral and (b) longitudinal correlation coefficients for 
$Z_{ref}=0$. Data points for the 0WT cases in panel (a) are mirrored around 
$y/D=0$. Grey lines represent the dimensions of the rotor, while the dashed lines indicate 
$\rho =0$.
Next, the longitudinal correlation coefficients are presented in figure 7(b). Here, the upstream probe is at the position indicated in the legend and 
$\Delta x/D$ indicates the non-dimensional separation distance to the second probe. This figure shows that the presence of the turbines disrupts the longitudinal coherence in the boundary layer indicated by a significantly lower correlation coefficient. Even at 
$x=11.8D$ downstream of the single turbine, the correlation is yet to recover to its unperturbed value. This reduction (of the order of nearly 50 % at hub height for 
$x=5.8D$) has been previously documented (Lukassen, Stevens & Wilczek Reference Lukassen, Stevens and Wilczek2017). The presence of the second turbine acts to further lower the correlation and increase the wake recovery length, as can be seen by data for 
$x/D=7$. The strength of the reduction of the longitudinal correlation in the wake flow is also affected by the reference height (omitted here), as previously discussed for the undisturbed flow in § 3.1.2.
In contrast to the observations in § 3.1.2, where the reduction of the turbulence length scales due to the thermal stability was monotonic and inhomogeneous in 
$L_x, L_y, L_z$, the turbines’ wakes have a comparable effect in all directions, where reductions in the length scales are of comparable size. Furthermore, the siting of the turbine within the array also influences the shape, and magnitude, of the correlation coefficients.
3.2.2. Effect of thermal stability on mean profiles and turbulence statistics
Vertical profiles in the wake of single- and two-turbine cases are presented in figures 8(a) and 8(b), respectively. Figure 8(a) shows that the far-field longitudinal wake recovery is impeded by increased thermal stability (see solid, dashed and dotted lines for 
$x/D=11.8$) as seen by the larger velocity deficit for the two stable cases; this is linked with the reduced turbulent activity within the incoming wind, and the reduction in mixing (Hancock & Pascheke Reference Hancock and Pascheke2014b; Machefaux et al. Reference Machefaux, Larsen, Koblitz, Troldborg, Kelly, Chougule, Hansen and Rodrigo2016). In contrast, the near-wake region is largely unmodified by the stability conditions (see profiles for 
$x/D\leqslant 2$). This is in line with findings highlighted by Hancock & Pascheke (Reference Hancock and Pascheke2014b) who noted that the ‘early part of the wake…develops independently of stratification’. The largest differences are found below hub height; this is interesting as it was shown how a change of inversion conditions leaves the bottom third of the boundary layer largely unaffected (Hancock & Hayden Reference Hancock and Hayden2020). The two-turbines case is presented in figure 8(b). A noticeable difference is the apparent more significant effect of thermal stability particularly in the near wake (compare red profiles at 
$x/D=2$ in figure 8a with blue counterpart at 
$x/D=8$ in figure 8b). However, this is due to the more persistent wake of the upstream turbine impinging on the second turbine in stable cases, consistent with the notion that for stable conditions, the turbulence transport is decreased in the wake flow (Du et al. Reference Du, Ge, Zeng, Cui and Liu2021), with a slower momentum recovery.
Figure 8. Vertical velocity component 
$U/U_e$ for (a) 1WT and (b) 2WTs. Data for 0WT are also included in black. Grey lines represent the dimensions of the rotor.
Complementary second-order velocity profiles (to figure 8) are shown in figure 9. For clarity, we only consider locations at 
$x=5.8D$. Reference cases (undisturbed boundary layers) in all three stability conditions are now included to aid the discussion. Well-known wake features are visible (see further discussion of figure 14). Consistent with the concept of ‘active filtering’ of Chamorro et al. (Reference Chamorro, Guala, Arndt and Sotiropoulos2012) is the enhancement of turbulence fluctuations in the upper part of the wake (
$Z/D>0.2$ in figure 9a) which is due to the blades rotation and associated tip vortical structures and their reduction (when compared to the undisturbed flow) in the lower part of the wake due to the disruption of the large-scale coherence of the boundary layer (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009), and the changes imposed to the mean shear. This is visible regardless of the thermal conditions of the boundary layer (compare red/blue dashed and dotted lines with the respective baseline cases in black). The main effect of stable stratification is the attenuation of the turbulence across the entire wall-normal range (Abkar & Porté-Agel Reference Abkar and Porté-Agel2015; Du et al. Reference Du, Ge, Zeng, Cui and Liu2021). The stability effects are much stronger in the top part of the wake (
$Z/D>0$) than in its lower counterpart. It was just discussed how the wake flow enhances the turbulence levels in the top of the boundary layer while attenuating them closer to the wall, then the thermal stability of the boundary layer acts in the same and opposite manner in the bottom and top part of the boundary layer, respectively.
Figure 9. Vertical velocity profiles 
$\overline {u^2}/{U_e}^2$ for (a) 1WT and (b) 2WTs. Data for 0WT are also reported in black. Grey lines represent the dimensions of the rotor.
3.2.3. Effect of thermal stability on velocity correlations
Having discussed the effect of wind stability on single-point measurements in the undisturbed flow (§ 3.1.1) and that of wind turbines on two-point measurements in neutral flow (§ 3.2.1) separately, this section investigates their effects in combination. Figure 10 shows the vertical coherence for the case of a single turbine, where both the presence of the turbine and the thermal stratification have a detrimental effect on the extent of the vertical coherence. The bottom part of the boundary layer (
$Z/D<-0.35$ in figure 10a) is more affected by the wake decay processes than by the thermal stability. The reduction in correlation coefficient is 12 % for a change in stability conditions (stable to mid-inversion), while approximately 50 % in the wake flows for neutral cases. This trend is similar at hub height in figure 10(b). However, the situation is different in the top layer (figure 10c for 
$Z/D>0.35$), where changes are more significant due to the stability conditions than to the presence of the turbine alone. Similar plots to those in figure 10 are shown for the wake-impinged turbine in figure 11. Here, changes are more visible in the near wall (figure 11a) rather than further up in the boundary layer (figure 11b,c.) Strikingly, almost no effect of stability is seen in figure 11(c). The temperature inversion has very mild effects on the vertical length scales close to the surface (figure 11a), where the dashed and dotted lines are almost indistinguishable, while more significant variations are found to be present at hub height in figure 11(b). Here, that the introduction of the second wind turbine dominates the flow physics rather than the strength of the thermal stability is visible by comparing profiles upstream (symbols *) and downstream (symbols 
$\bigtriangledown$) of the turbine. Interestingly, the presence of the second turbine seems to facilitate the recovery of the wall-normal coefficients above hub height (
$Z/D>0$) – perhaps fortuitously.
Figure 10. Vertical correlation coefficients in the presence for 1WT for (a) 
$Z_{ref}=-0.35D$, (b) 
$Z_{ref}=0$ and (c) 
$Z_{ref}=0.35D$. Data for 0WT are also reported in black. Grey lines represent the dimensions of the rotor, while the vertical dashed lines indicate 
$\rho =0$.
Figure 11. Vertical correlation coefficients in the presence for 2WT for (a) 
$Z_{ref}=-0.35D$, (b) 
$Z_{ref}=0$ and (c) 
$Z_{ref}=0.35D$. Data for 0WT are also reported in black. Grey lines represent the dimensions of the rotor, while the vertical dashed lines indicate 
$\rho =0$.
Finally, the correlations in the other two directions are presented in figure 12. In figure 12(a), the wake flows reduce the lateral correlation whilst also enhancing its asymmetry (as previously discussed in § 3.2.1); this is particularly true for the case with multiple wind turbines (see blue lines compared with black). The wind stability is acting similarly to what has been described so far, where stabler winds tend to reduce the length scales; its effect is more pronounced in the vertical direction than it is on the lateral coherence. Also note that, overall, the thermal stability opposes the asymmetry in the correlation caused by the turbines (see the reduction of negative tails for both positive and negative 
$y$), as suggested by Hancock & Pascheke (Reference Hancock and Pascheke2014b). In this respect, the presence of the turbines and the thermal stability have competing effects. It is also interesting to note the behaviour of the longitudinal correlation coefficient in figure 12(b). Here, the presence of the turbines and the thermal stability both contribute towards a significant shortening of 
$L_x$. We ascribe this, in the former, to the disruption of the coherent boundary layer structures due to the wake flows, while in the latter, to the suppression of vertical momentum. This reduction in the longitudinal coherence is much more sensitive to the presence of the turbines than to the thermal stability.
Figure 12. (a) Lateral and (b) longitudinal correlation coefficients in different thermal conditions. Data points for the 0WT cases in panel (a) are mirrored around 
$y/D=0$. Measurements in panel (b) are taken at hub height. Grey lines represent the dimensions of the rotor, while the dashed lines indicate the 
$\rho =0$ location.
3.2.4. Effect of thermal stability on spectra
Finally, this section discusses the effect of thermal stability on the power spectral density within the near-wake flow. Results are presented in figure 13. The most notable feature in the near field is the distinct signatures of the rotor and blade-passing frequency (in the form of sharp spikes in the spectra) visible in correspondence of the top edge of the rotor in figure 13(b) downstream of the first turbine for the stable cases (broken red lines). The rotor frequency 
$f_{wt}$ is expected to be in the range 
$5 {\rm Hz} < f_{wt}<6.6\ {\rm Hz}$ with the blade-passing frequency in the range 
$15\ {\rm Hz} < f_{bpf}< 20\ {\rm Hz}$ (highlighted by the shaded area). The spectra for the stable cases feature peaks which correspond well with these frequencies (particularly considering the uncertainty in the rpm controller). Energy peaks are also visible at harmonics and sub-harmonics of the rotational frequency; these have been previously associated with vibrations, rotor mass imbalance and mechanical looseness at the base of the turbine (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2017). These features are absent in the neutral case – an indication of the much more persistent signature of the blade-induced turbulence in stable thermal conditions. It should be noted that the tip vortices leave a stronger signature in the vertical velocity component (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2010), which therefore would be interesting to measure. The blade-passing frequency signature is not identifiable downstream of the second turbine, where the higher background turbulence levels help dissipate the tip vortices much closer behind the rotor. Lastly, figure 13(a) strengthens the view of the turbine as an ‘active filter’ (see Chamorro et al. (Reference Chamorro, Guala, Arndt and Sotiropoulos2012) and Appendixes A and B), while the thermal stability indiscriminately lowers the energy at all turbulent scales. The effect of the temperature inversion is minimal at this wall-normal height for both turbine cases. As for the far wake (omitted here), findings are similar to those described in figures 15(b) and 15(d) in Appendix B. To summarise, a stability modulation of the spectral content of the turbulence is visible across the whole frequency range, but more significantly at low frequencies confirming that thermal effects involve predominantly large flow scales. Within the wake flows, the temperature inversion has a highly localised effect, with differences in energy which are discernible across all flow scales, but not as significant as in those recorded for the undisturbed case (figure 3e, f). This hints at the complex and somewhat juxtaposed nature of the interaction between thermal and wake effects. The latter energises the small scales and hinders the near-wall large-scale coherence of the turbulence; the former suppresses the turbulent fluctuations, while also rendering the wake flow more persistent.
Figure 13. Power spectral density of streamwise velocity fluctuations at 
$x/D=1$. (a) 
$Z=-0.5D$ and (b) 
$Z=0.5D$. All data are taken at 
$y/D=0$. The spectral slope represents the inertial subrange (
$\beta =-5/3$).
4. Conclusions
Wind tunnel tests were performed in both neutrally and moderately thermally stratified boundary layers (with and without an overlying inversion) to study the effects of the thermal stability on turbulent properties in the wake flow of wind turbines. A baseline case (undisturbed boundary layer) and cases with one and two turbines in line were considered, for a total of nine different scenarios (when including the different thermal stability states). Streamwise velocity and two-point velocity correlations were acquired via means of two independent LDA probes, which allow one to characterise mean and fluctuating velocities, and to undertake frequency analysis. The analysis emphasises the role of both the presence of the wind turbines (and their wake flows) and the thermal stability in significantly influencing the flow physics.
In the undisturbed boundary layer, the thermal stability affects the turbulence length scales predominantly in the vertical and longitudinal directions (by 11 % and 23 % in the middle of the layer compared to the neutral counterpart), while changes are largely negligible in the lateral direction. This is due to the selective dampening of the vertical momentum transport in stable boundary layers (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017; Hancock & Hayden Reference Hancock and Hayden2018, Reference Hancock and Hayden2020). The introduction of a temperature inversion leaves the turbulence length scales in the bottom third of the boundary layer largely unaffected while the top of the layer sees further reductions (to 51 % and 43 %). This also translates to the spectral content of the turbulence: in the near-wall region, it is indistinguishable with and without the overlying inversion, while further up in the boundary layer, the more significant differences are those affecting the large scales of the flow. These findings extend the observation of Hancock & Hayden (Reference Hancock and Hayden2020) to velocity correlations and spectral content, both of which were not considered before.
In contrast to what was observed for thermal effects, the introduction of wind turbines (and their associated wake flows) in neutral boundary layers results in a similar reduction of the length scales in all three directions (50 %, 40 % and 49 % at 
$x/D=5.8$ and hub height) compared to the undisturbed case, with asymmetries in the lateral correlation arising due to the rotor-induced wake rotation. The second turbine helps the vertical and lateral correlation recovery (particularly in the outer part of the boundary layer), while contributing to further reducing its longitudinal counterpart.
The effect of different thermal stabilities on the turbulence length scales in the wake flow is complex. In the bottom part of the boundary layer, the vertical correlation is much more affected by the wake flow (reduced by 
$\approx$50 %) than by the changes in thermal stability (where the reduction is of the order of 
$12\,\%$). The situation is reversed in the top half of the layer, where the thermal stability is the main culprit for the length scale reduction. Effects are also evident on the longitudinal coherence, which is significantly reduced as a result of both the wake flows and the thermal stability (the first factor being the dominant one). Finally, blade rotation is also known to introduce asymmetry in the lateral coherence within the wake turbulence, while the thermal stability is here shown to seemingly oppose this effect. As for thermal modulation of the turbulence characteristics in the wake flow, the stability acts to inhibit the turbulent fluctuations, while also simultaneously increasing the wake flow persistence, in accordance with most of the previous literature on this topic (Magnusson & Smedman Reference Magnusson and Smedman1994; Hancock & Pascheke Reference Hancock and Pascheke2014b), but not Chamorro & Porté-Agel (Reference Chamorro and Porté-Agel2010). Furthermore, this work confirms that the wake-impinged turbine typically shows a higher momentum deficit, though its wake recovery rate can be faster due to the higher turbulence levels present in the wake. Interestingly, given the known effect of wind turbines on turbulence levels (enhancement/attenuation in the top/bottom half of the wake), the thermal stability acts in the same and opposite manner in the bottom and top part of the boundary layer, respectively.
Although the cases of only one and two turbines in line have been considered herein, findings suggest that, for larger arrays, the siting of the turbines within the array (and presumably within the farm) also can influence the shape and magnitude of the correlation coefficients and their recovery. This would indicate different fatigue loading across the farm. Furthermore, our results strongly suggest that both the effect of thermal stability and the local flow characteristics due to turbine siting should be considered – in combination – when considering the application of Veers-type approaches to stable winds, as their interplay contributes to modify the turbulent correlation maps, and hence the wind coherence functions in all three directions. Further work should be expanded to larger turbine arrays to confirm these findings.
Acknowledgements
The authors want to acknowledge Mr A. Makedonas for his initial assistance with the data organisation.
Funding
The work herein was done with support from the EPSRC under the MAXFARM project (EP/N006224/1) and via EPSRC-Flexible Funding agreements within SUPERGEN Wind (coordinated by the University of Strathclyde) and Supergen ORE (coordinated by the University of Plymouth) under projects VESABL and VENTI, respectively. The authors are also grateful to NERC/NCAS(AMOF) and EPSRC for supporting the EnFlo wind tunnel as a national facility.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in Open Research at http://doi.org/10.15126/surreydata.900776.
Author contributions
The experiment was designed by M. Placidi and P.E. Hancock. The tests were led by P. Hayden, who also coordinated the data collection, and created and administered the data files. Data reduction, figure preparation and the first draft of the manuscript were carried out by M. Placidi. All authors reviewed the manuscript.
Appendix A. Mean and fluctuating velocity profiles in neutral wake flows
The streamwise wake development in the case of single and two wind turbines is documented in figures 14(a,b) and 14(c,d), respectively. It can be seen in figures 14(a,c) that, in the near wake (
$1D< x<5.8D$), the profiles display the classical two-lobe shape, where areas of momentum deficit are visible over the blade length (Hancock & Pascheke Reference Hancock and Pascheke2014b). The wake does not fully recover by 
$x=11.8D$, where the velocity profiles in the wake flow do not match that of the undisturbed flow (black line). The presence of the turbines also visibly influences the flow upstream – compare the red *-profile at 
$x=-0.2D$ with the undisturbed flow (black line) in figure 14(a) and the blue *-profile at 
$x=5.8D$ with the red 
$\triangle$-profile in figure 14(c). This is due to the deceleration associated with the impedance of the turbines. The wake of the second turbine in figure 14(c) behaves similarly, however, some of the features across the rotor are less pronounced and a single-lobe wake establishes earlier, due to its waked position. This wake also shows signs of stronger persistence where the momentum deficit compared to the undisturbed conditions is much larger for 
$x=11.8D$ in figure 14(c) than for 
$x=5.8D$ in figure 14(a). The presence of the second turbine is also visible in the profiles just upstream of it (
$x=5.8D$) in figure 14(c). The lateral profiles (not shown) ahead of the first turbine are found to be symmetrical, while those ahead of the second turbines are not; this can be ascribed to the wake flow rotation.
Figure 14. Vertical velocity profiles (a,c) 
$U/U_e$ and (b,d) 
$\overline {u^2}/{U_e}^2$ for (a,b) 1WT and (c,d) 2WTs. Data for 0WT are also reported in black. Grey lines represent the dimensions of the rotor.
The streamwise velocity fluctuations are presented in figures 14(b) and 14(d). Only a limited number of streamwise stations are discussed, as these are used to exemplify the most salient flow features. First, the profiles are shown to be asymmetric with respect to the hub location; this asymmetry is due, of course, to the mean vertical shear and the rotation of the blades (particularly obvious from lateral profiles omitted here). As previously discussed, the enhancement/attenuation of turbulence fluctuations in the upper/lower part of the wake is consistent with previous findings (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009; Chamorro et al. Reference Chamorro, Guala, Arndt and Sotiropoulos2012). The imprint of the wake is still visible in the turbulence fluctuations at 
$x/D=11.8$ in figure 14(b), although less pronounced than in the near field at 
$x/D=1$. Similar remarks can be made for the wake-impinged turbine (figure 14d); however, the enhancement of the turbulent fluctuations is even more pronounced here, due to the wake-flow superposition. One further interesting difference between the upstream and wake-impinged turbines is within their respective induction zone (i.e. their upstream effect). In the upstream turbine, turbulent fluctuations in the induction region are enhanced across the entire depth of the boundary layer (compare the black line with the red *-profile in figure 14(b); this is in line with the findings of Li et al. (Reference Li, Huang, Ge and Zhang2022). In the wake-impinged turbine, instead, the turbulent fluctuations in the induction zone are reduced across the whole wall-normal range (compare the red 
$\triangle$-profile with the blue *-profile in figure 14d). This opposite effect is presumably due to the different flow characteristics (and length scales) to which the wake-impinged turbine is subject. Overall, when comparing the fluctuations upstream of the second turbine with those in the undisturbed boundary layer, these are still augmented by the wake-added turbulence for 
$Z/D>-0.2$, but the opposite is true below this point. This is a combination of the upstream turbine's active filtering effect, but also a genuine – though perhaps counterintuitive – induction effect of the second turbine.
Appendix B. Spectra in neutral wake flows
Figure 15 displays the effect of the presence of wind turbines on the power spectral density of the streamwise velocity fluctuations. Both data from the near and far fields (
$x/D=1$ and 
$x/D\approx 6$) are presented in figures 15(a,b) and 15(c,d), respectively, while two different wall-normal heights are considered (corresponding to the edge of the rotor). In the near-field data in figure 15(a), the effect of the blade rotation and associated tip vortical features are identifiable by higher energy at the high frequencies, hence the smallest flow scales. This is due to the turbine's activity as an active turbulence generator (Chamorro et al. Reference Chamorro, Guala, Arndt and Sotiropoulos2012). Naturally, changes are more noticeable at wall-normal locations corresponding to the top edge of the rotor (figure 15c) than at the bottom in figure 15(a) (due to a combination of the mean shear and the presence of the turbine tower Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2009, Reference Chamorro and Porté-Agel2010); these features are stronger downstream of the second turbine. It should be noted that for these cases, the blade-passing frequency is between 15 and 18 Hz, and no strong peak is visible at this frequency for 
$x/D=1$; though this is expected (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2010; Chamorro et al. Reference Chamorro, Arndt and Sotiropoulos2011; Hu, Zifeng & Partha Reference Hu, Zifeng and Partha2012). There is also no discernible spectral feature of the shedding induced by the rotor. The lack of the blades’ rotation signature in the spectra in figure 15(b) is believed to be due to the relatively high level of turbulence intensity, 
$Tu=1/U_\infty \sqrt {1/3(\overline {uu}+\overline {vv}+\overline {ww})}=4.23\,\%$, characterising our wind tunnel, which quickly acts to dissipates the blades’ tip vortices. This was previously suggested by Hancock & Pascheke (Reference Hancock and Pascheke2014b) and it is further explored in § 3.2.4. The turbulence levels here are considerably higher than customary values reached in facilities with flow conditioning, though still relatively low compared with much of the previous experimental work on this topic (e.g. Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014). Also of note in figure 15(a) is the attenuation of the energy at low frequency in the wake flow; this has previously been described as the turbine acting as an ‘active filter’ energising the high-frequency range while dampening the low frequencies (Chamorro et al. Reference Chamorro, Guala, Arndt and Sotiropoulos2012) closer to the wall. This feature is absent for 
$Z=0.5D$ in figure 15(c). Figures 15(b) and 15(d) show correspondent data to those just discussed but at 
$5.8D$ downstream of the turbines. Findings are as expected, with the upper wake flows showing higher energy densities as the number of wind turbine increases and their wake flows superimpose for all wavenumbers. The energy contained in the fluctuations for the undisturbed case (black line) is higher for 
$Z=-0.5D$ in figure 15(b), in accordance with turbulence wake profiles in figures 14(b and 14(d), which appear larger in this case for 
$Z/D<-0.3$ compared to the case with the turbines. Although not shown here, the hub-height spectrum at 
$x/D=11.8$ downstream of the first turbine is far from recovered to the undisturbed case.
Figure 15. Power spectral density of streamwise velocity fluctuations at (a,c) 
$x/D=1$ and (b,d) 
$x/D\approx 6$. (a,b) 
$Z=-0.5D$ and (c,d) 
$Z=0.5D$. All data are taken at 
$y/D=0$. The spectral slope represents the inertial subrange (
$\beta =-5/3$).
Appendix C. Effect of thermal stability on wake recovery
To further investigate the wake recovery rate, we apply the framework of Magnusson & Smedman (Reference Magnusson and Smedman1994) and Barthelmie et al. (Reference Barthelmie2004). The local momentum deficit, 
$\Delta U$, is normalised by the free stream wind velocity at hub height, 
$U_{h}$, and then expressed as a power-law behaviour as
\begin{equation} \frac{\Delta U}{U_{h}}=K\left(\frac{x}{D}\right)^n, \end{equation}
where 
$K$ is a constant (function of the atmospheric stability only) and 
$n$ is nominally invariant – herein, 
$n=-0.8$. The momentum deficit is calculated from lateral horizontal profiles at hub height measured across the whole wake. Results for both the single and tandem turbine cases are presented in figures 16(a) and 16(b), respectively. Here, best-fitting yielded 
$K=1.37$ (neutral) and 
$K=1.48$ (stable) in figure 16(a), while 
$K=1.32$ and 
$K=1.35$ (neutral and stable) in figure 16(b). These findings are consistent with observations in figure 8, where the wake recovery rate was found to be slower with increasing thermal stability. The momentum deficit is higher for the wake-impinged turbine (consistent with a superposition of effects), however, its recovery rate is faster (consistent with a wake flow characterised by higher turbulence levels). For comparison, models suggested by Magnusson & Smedman (Reference Magnusson and Smedman1994), Hancock & Pascheke (Reference Hancock and Pascheke2014b) and Barthelmie et al. (Reference Barthelmie2004) are reported as a bold black line and dark- and light-shaded areas, respectively. These alone already show the spread in the data from previous work. The best-performing fit is that of Hancock & Pascheke (Reference Hancock and Pascheke2014b), however, the agreement is far from perfect. Furthermore, these findings indicate that for the upstream turbine, thermal effects are perhaps more significant, while the wake-impinged turbine response is largely dominated by the wake physics. It is also interesting to note from figure 16(a) that in the near wake (
$x/D<2$), the rate of increase of momentum deficit seems to be more gradual for stable cases than in the neutral case. This is consistent with observations by Hancock & Pascheke (Reference Hancock and Pascheke2014b) who noted that atmospheric stability ‘…would give a reduced rate of blade wakes merging and mixing’, hence influencing the early wake development. Some caution is needed when interpreting these results due to both experimental and fitting uncertainties, however, this information is largely absent in the previous work (Abkar & Porté-Agel Reference Abkar and Porté-Agel2015; Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2010; Radünz et al. Reference Radünz, Sakagami, Haas, Petry, Passos, Miqueletti and Dias2021). Overall, results are qualitatively in agreement with most previous laboratory and field experiments which found a slower wake recovery in the far field due to thermal stability (Magnusson & Smedman Reference Magnusson and Smedman1994; Hancock & Pascheke Reference Hancock and Pascheke2014b). In contrast, Chamorro & Porté-Agel (Reference Chamorro and Porté-Agel2010) suggested that the velocity recovery appeared faster in stable boundary layers for 
$x/D<10$, while largely independent of thermal effects beyond this location – it is unclear why.
Figure 16. Maximum lateral velocity deficit normalised by 
$K$ as a function of non-dimensional distance from the turbines. (a) 1WT case and (b) 2WT case. Light- and dark-grey areas represent fitting from Barthelmie et al. (Reference Barthelmie2004) and Hancock & Pascheke (Reference Hancock and Pascheke2014b), while the black bold line indicate fitting by Magnusson & Smedman (Reference Magnusson and Smedman1994). Here, 
$x_r$ in panel (b) indicates the relative distance from the second turbine. Data from all three thermal states are reported.
Open access
