Impact of Negative Geostrophic Wind Shear on Wind Farm Performance

Baroclinicity, which leads to height-dependent driving pressure gradients, occurs in various situations such as the ﬂow transition between land and sea, and sloping terrain. It has been shown that baroclinicity modiﬁes the structure of the atmospheric boundary layer. For example, negative shear baroclinicity creates additional turbulence at higher elevations, which might inﬂuence the energy entrainment into large wind farms. Here, we use large-eddy simulations to study the eﬀect of baroclinicity-induced negative shear on the wind farm power production and energy entrainment into a large wind farm. In agreement with the literature, our simulations show that negative geostrophic wind shear signiﬁcantly modiﬁes the mean wind velocity in the atmospheric boundary layer. Speciﬁcally, for the cases considered in the study, the negative geostrophic shear causes a change in the mean velocity up to 2.3 m/s at hub height, which greatly alters the wind farm power production. Additionally, we demonstrate with an energy budget analysis that a wind farm does not necessarily beneﬁt from the additional turbulence created by the negative geostrophic wind shear. The reason for this is that the baroclinicity-induced negative shear alters the height and strength of the low-level jet and creates an upward ﬂux above the jet, limiting the energy entrainment into the wind farm. Our results show that wind resources are altered in the boundary layer due to negative geostrophic wind shear and should be considered in wind farm modeling and power forecasts.


I. INTRODUCTION
One of the uncertainties in wind farm power production is its dependency on the prevailing atmospheric conditions. For instance, energy entrainment into the turbine wakes from the atmosphere above plays a dominant role in the overall efficiency of a wind farm [1,2]. Consequently, understanding the interaction between the atmospheric boundary layer (ABL) and wind farms is instrumental in increasing the wind farm's efficiency. However, the variety of phenomena occurring in the ABL makes its description and modeling highly complex. Some of these phenomena are variations in atmospheric stability [3][4][5][6][7][8][9], cloud formations [7,10], and geostrophic wind [11][12][13][14][15][16][17][18]. In this study, we focus on the latter.
According to the thermal wind balance [19], horizontal temperature gradients can cause a variation in the geostrophic wind with height, which is known as baroclinicity [20]. A height-dependent geostrophic wind implies the friction velocity, Obukhov length, shear production, and ABL height strongly depend on baroclinicity when weakly stable conditions prevail. For stronger atmospheric stability, when the lower part of the ABL is decoupled from higher elevations due to the strong turbulence destruction, the dependency of these parameters on baroclinicity decreases [27]. Previous experiments and simulations have shown dependencies of the wind profiles on baroclinicity for neutral and unstable conditions as well [18,23,24,[28][29][30][31][32][33]. Recently, Ghannam and Bou-Zeid [34] introduced a correction to the logarithmic law predictions to account for baroclinicity effects in near neutral conditions. They used first-order closure principles to capture wind turning and baroclinicity effects.
As baroclinicity affects the structure of the ABL, it is expected to have a large impact on the performance of wind farms. The height and strength of the LLJ, which are highly influenced by baroclinicity, play an important role in the wake recovery inside wind farms [35][36][37]. Conangla and Cuxart [20] measured and simulated cases of negative geostrophic wind shear, where the LLJ is positioned at heights starting from 50 m and the turbulence kinetic energy increased above, i.e., between 75 and 250 m. Wind turbine rotors are typically positioned within this range of heights. The increased turbulence intensity that results from negative geostrophic wind shear is expected to influence the momentum transport in the ABL, and hence the ability of the wind farm to harvest energy from it.
However, how geostrophic shear influences the entrainment into wind farms is still unknown. To study the complex interaction of wind farms and turbulent ABL flows, Calaf et al. [38] proposed using large-eddy simulations (LESs). The advantage of LESs is that they capture temporal fluctuations and resolve large-scale flow features in the ABL, while the small-scale turbulence is parameterized using a subgrid scale (SGS) model [22]. Thus, LESs can accurately model the complex flows through wind farms. The first LESs of wind farms assume neutral pressuredriven flow over flat terrain [38] and this simplification is still commonly used in the fundamental research on wind farm flows [39][40][41][42][43][44]. Studying flows under these controlled conditions is necessary to get more fundamental insight. While LESs of wind farms in neutral pressure-driven flow have been shown to reproduce measurements [41], their results are not universally applicable. The reason for this is the underlying assumption in such simulations that the wind turbines reside in the inner region of the ABL and that outer layer effects are negligible. In later studies, outer layer effects like thermal stability [45][46][47][48][49][50][51][52][53] and Coriolis force [54][55][56] received more attention. These studies highlighted the complexity of the interaction between wind farms and the ABL. In particular, they showed that the wake recovery and the energy entrainment into the wind farms are highly dependent on the atmospheric conditions, such as the height and strength of the capping inversion and the free-atmosphere stratification. Because of the high complexity of the interaction between ABL dynamics and wind farms, understanding the influence of atmospheric mesoscales on microscale processes was identified as one of the grand challenges in wind energy [57].
Recently, realistic flow conditions have been included in LESs via data assimilation [58,59], or by coupling LESs to mesoscale models like the Weather Research and Forecasting model [60,61]. However, the complexity involved in these approaches makes it difficult to isolate and understand the effect of baroclinicity on wind farm performance. The present study aims to identify the differences in wind farm power production with and without a prevailing negative geostrophic shear. We especially want to find out how the altered velocity profile and the increased turbulence intensity due to negative shear baroclinicity affect the energy entrainment into the farm; see Fig. 1.
The remainder of this article is structured as follows. Section II discusses the simulation approach and the considered cases. The effect of negative geostrophic wind shear on atmospheric dynamics is addressed in Sec. III. The effect of negative geostrophic wind shear on the flow in and around wind farms is studied in Sec. IV, followed by an analysis of the entrainment fluxes (Sec. V) and the effect on the power production (Sec. VI). We end with the conclusions and outlook in Sec. VII.

II. LARGE-EDDY SIMULATIONS
We perform LESs with an updated version of the code used by Albertson and Parlange [62]. The code was validated by Gadde et al. [63] and can accurately simulate thermally stratified ABLs and wind farm wakes [64]. The governing equations for the LESs, in terms of the filtered quantities (denoted by a tilde), are the incompressible continuity equation, the equation of momentum conservation, and the transport equation for potential temperature: Here, i = 1, 2, 3 correspond to the streamwise (x, u), spanwise (y, v), and vertical (z, w) directions, respectively. Furthermore,ũ i represents the filtered velocity field components andθ is the filtered potential temperature field. The Boussinesq approximation is applied to obtain the buoyancy term in Eq. (2), with gravitational constant g, Kronecker delta δ ij , and buoyancy parameter β = 1/ θ with 023007-2 respect to the planar-averaged potential temperature θ . Effects of resolved viscous stresses are neglected, since a very high Reynolds number flow is assumed. The SGS stress tensor is denoted by τ t ij and τ ij = u i u j −ũ iũj is its deviatoric, traceless part, while the trace of the SGS stress tensor τ t kk /3 is absorbed into the filtered modified pressurẽ p * =p/ρ 0 − p ∞ /ρ 0 + τ t kk /3 with air density ρ 0 . The SGS heat flux vector is given by q i = u i θ −ũ iθ . Both SGS deviatoric stress and SGS heat flux are modeled using an eddy diffusivity parameterization. Specifically, the anisotropic minimum dissipation model [63,65] is applied. The model provides dynamic and scale-dependent SGS coefficients for modeling SGS turbulence. The geostrophic wind velocity is given by G i = − ij 3 ∂ j p ∞ /(ρ 0 f c ), where ij 3 denotes the alternating unit tensor; f c is the Coriolis parameter; and ∂ i p ∞ is the driving mean pressure gradient, which is independent of height in barotropic conditions, but depends on height in baroclinic conditions. The computational domain is discretized with n x , n y , and n z points in the streamwise, spanwise, and vertical directions, respectively. A uniform grid is used in horizontal directions, with a corresponding grid spacing of x = L x /n x and y = L y /n y , where L x and L y are the dimensions of the computational domain. In the vertical direction the grid is uniform (with spacing z ) up to a height z uni . Beyond z uni the grid is stretched up to the domain height L z , using a hyperbolic tangent stretching function. The computational grid is staggered in the vertical direction. The variables u, v, and θ are stored at the intervals z stag l = (z l + z l+1 )/2, with l = 0, . . . , n z − 1, whereas the vertical velocity w is stored in the nodes z l .
Free-slip boundary conditions with zero vertical velocity are enforced at the top of the domain. A Rayleigh damping layer is used at the top of the domain to reduce the reflection of gravity waves that are triggered by the presence of the wind farm [66]. For the bottom boundary condition, the wall shear stress τ i3|w and the buoyancy flux q * at the surface are modeled using the Monin-Obukhov similarity theory [5]: Here u * is the frictional velocity, z 0 is the roughness height, κ = 0.4 is the von Kármán constant,ũ r = √ũ 2 +ṽ 2 is the filtered velocity magnitude at the first grid level [67], θ s is the filtered potential temperature at the surface, and z 0,s = z 0 /10 [68] is the thermal surface roughness length. For stable boundary layers (SBLs), the stability corrections for momentum ψ M = −4.8z/L and heat flux ψ H = −7.8z/L are used [69]. Here θ 0 is the reference potential temperature and L = −(u 3 * θ 0 )/(κgq * ) is the surface Obukhov length. For neutral boundary layers (NBLs), the buoyancy flux at the surface and the stability corrections are zero (q * = 0 Kms −1 , ψ M = ψ H = 0). Time integration is performed using a third-order accurate Adams-Bashforth scheme. Spatial derivatives in the vertical direction are calculated using a second-order central finite difference scheme. A pseudospectral method is applied in horizontal directions, resulting in periodic boundary conditions in horizontal directions. The concurrent precursor method [70] is employed to generate realistic, turbulent inflow conditions. This approach samples flow data from a periodic ABL simulation performed in a precursor domain. The sampled data are introduced as inflow velocity into a fringe region of the wind farm domain using a symmetric fringe function [71]. As the wind direction changes with height, a second fringe region is employed in the spanwise direction. We employ a proportional integral controller [72] to guarantee that the planar-averaged wind angle at hub height is 0 • . This is necessary as the wind veer depends on the stability conditions.

A. Wind turbine parameterization
The turbine force term f i is implemented using the actuator line method [73,74]. The turbine blades are represented by distributed body forces, calculated dynamically using the local flow velocity. The total force f i comprises a lift and drag component, evaluated in a local coordinate frame (r, θ, x) before being transformed to the global frame (x, y, z). The local lift and drag forces, per unit span, are given by respectively. Here,ũ rel = ũ 2 θ +ũ 2 x is the flow velocity relative to the blade, with circumferential velocityũ θ = r −ũ y cos(θ ) +ũ z sin(θ ). Furthermore, is the rotational speed of the rotor and c is the local chord length. The lift and drag coefficients C L (α) and C D (α), at the local angle of attack α = φ − γ , are obtained from tabulated airfoil data. Here, φ is the angle betweenũ θ andũ rel , while γ accounts for any twist and pitch contributions. A Gaussian projection, given by η = −3 π −3/2 exp(−d 2 / 2 ), where denotes the kernel width and d is the distance to the considered actuator point, is used to smear the turbine force and avoid numerical instabilities. The kernel width for the force projection is related to the grid spacing via = 2 x , based on initial testing, as well as the recommendation by Martínez-Tossas et al. [75]. For practical reasons, the rotor heads are oriented perpendicular to the x direction rather than perpendicular to the local flow. In the Appendix, we show that this does not influence the main findings of the study.
We note that the LES filtering tilde is omitted to simplify the notation in the remainder of this work.

B. Modeling baroclinicity in LESs
To evaluate the effects of baroclinicity on wind farm performance and wake development, we study cases with varying geostrophic wind profiles and boundary-layer stability. Baroclinicity is generated by temperature gradients, resulting in thermal winds within the ABL, which weaken aloft [24]. While the temperature gradients occur on very large horizontal scales, the geostrophic wind changes with height also on small horizontal scales. For instance, a temperature change of 1 K over 1000 km in the horizontal direction can cause an increase in geostrophic velocity of 3 ms −1 km −1 with height [28]. On the relatively small domain size considered here, we model the effect of baroclinicity by varying the geostrophic wind with height. As pointed out by Sorbjan [33], this is a relatively straightforward approach.
The boundary layer is driven by a geostrophic wind with magnitude G. Cases with constant G are known as barotropic, while cases where G is height dependent are named baroclinic. The selection of the geostrophic wind velocity for the SBL cases in this study is inspired by the measurements of Conangla and Cuxart [20], who observed and modeled negative shear baroclinicity in a SBL. Interestingly, they observed that the negative geostrophic wind shear gives rise to sharp changes in the LLJs and increased turbulence intensity aloft. The changes in the shape of the LLJ and the increased turbulence intensity above the LLJ are expected to change the energy transport in the ABL. The entrainment from above highly influences wind farm performance. Consequently, we are interested in how the increased turbulence intensity observed by Conangla and Cuxart [20] affects this entrainment process. They used a turbulent kinetic energy single-column model to show that the measured wind profiles could result from the following baroclinic geostrophic wind profile: Here, G 0 refers to the geostrophic velocity at the surface. All baroclinic cases under consideration have a negative shear layer in the geostrophic velocity profile. It is modeled by a linear vertical variation in the geostrophic wind of magnitude G over a height of z s , starting from a height z s . In all cases, the geostrophic velocity is constant above an altitude of z s + z s ; see also Fig. 2 below. Additionally, we are interested in how negative geostrophic wind shear effects are different in SBLs and NBLs. Therefore, we draw inspiration from the study by Floors et al. [23] and simulations by Momen et al. [24] to design realistic NBL conditions with negative geostrophic shear, allowing us to study wind farm performance under these conditions. The NBL cases are also constructed using Eq. (7).

C. Suite of LESs
We consider a domain size of L x = 15.36 km, L y = 4.8 km, and L z = 4.0 km. Previous studies [64] confirmed that this domain size allows the turbulence statistics relevant to this study to be captured accurately. The domain is discretized using 1280 × 640 × 384 grid points in the streamwise, spanwise, and vertical directions, respectively. The streamwise and spanwise grid resolutions are x = 12 m and y = 7.5 m, respectively. The vertical grid spacing equals z = 5 m up to z uni = 1.5 km, and is thereafter stretched to a maximum of z = 59 m at the top boundary. The fringe region covers the final 1.54 km in the streamwise direction and 0.48 km in the spanwise direction. The roughness length is z 0 = 0.002 m, corresponding to offshore conditions [76] and the reference potential temperature is θ 0 = 286 K. The initial potential temperature profile contains a mixed layer with a constant potential temperature of θ 0 up to 1200 m. Starting from z = 1200 m, a capping inversion is implemented by increasing the potential temperature by 3K over a height of 200 m. The free atmosphere is positioned above z = 1400 m with a stratification strength of 5 K/km. Furthermore, we apply a constant surface cooling of C r = 0.25 K/h for the SBL cases, while the surface temperature is kept constant for the NBL cases. The Coriolis parameter is set to f c = 1.159 × 10 −4 s −1 (corresponding to latitude = 52 • ). The details of the geostrophic wind used in the study and its variation with height are listed in Table I and shown in Fig. 2, respectively. The magnitude of the geostrophic wind velocity for the barotropic cases is selected to achieve hub height velocities around 10 m/s, which are typical for offshore wind farms. The geostrophic wind velocities in the baroclinic cases match the barotropic velocities at the top or bottom. In this way, we can study how the actual velocity magnitude and LLJ height influence the result. The initial wind profile is set equal to G(z); see Eq. (7). Random perturbations are added to the velocity profile below 200 m to spin up turbulence. The amplitude of the perturbations has a maximum value of 10 −5 G 0 at the ground and decreases linearly with height. Similarly, random perturbations with a maximum magnitude of 10 −5 θ 0 are added to the initial potential temperature profile. The boundary layer is assumed to be in a quasisteady state when the velocity and other turbulent quantities have reached a steady state, and the temperature profile changes at a constant rate [77].
We perform the simulations in two stages. In the first stage, we perform a spin-up simulation in a domain of size L x /2 × L y /2 × L z for 7 h. In the second stage, after the flow has reached a quasisteady state, we use the domain periodicity to initialize the flow for the wind farm domain and activate the wind turbines. The statistics are collected over 3 h starting from the 8th hour to the 11th hour.
We consider a wind farm with 10 × 5 wind turbines. A schematic of the setup is shown in Fig. 3. The turbines are positioned uniformly in an aligned layout, with a spacing of s x = 7D and s y = 5D in the streamwise and spanwise directions, respectively. We simulate turbines based on the National Renewable Energy Laboratory (NREL) 5-megawatt turbine [78], which has a rotor diameter of D = 126 m and a hub height of z h = 90 m. The tip-speed ratio (TSR) is fixed to λ = 7.55, which gives an optimal power coefficient C p [78]. Consequently, the rotor rotational speed is computed dynamically to accommodate the prescribed TSR. Dependent on the local inflow conditions, varies from 6.7 to 13.4 rpm and the thrust coefficient lies between 0.84 and 0.87 in the simulations under consideration. Figure 4 shows the mean horizontal wind magnitude v g = u 2 + v 2 for the SBL (a) and NBL (e). Here, · and · denote planar and temporal averaging, respectively. The geostrophic forcing [Eq. (7)] is displayed with thin lines. The actual velocity matches the imposed geostrophic wind velocity G 0 − G at heights above the capping inversion (shaded area in Fig. 4). This is expected as above the BL height, the turbulent friction is small, and the geostrophic balance between the Coriolis and pressure gradient forces is reached. Below the capping inversion, the flow velocity strongly depends on the imposed geostrophic shear. Specifically, the jet strength depends on the local forcing strength; see also Table II. The jet height is lower for the baroclinic cases than for the barotropic cases. Similar observations were made by Momen et al. [24] and Conangla and Cuxart [20], who also simulated and measured geostrophic velocities with negative shear. To conclude, Figs. 4(a) and 4(e) show that negative shear baroclinicity greatly alters the mean wind profile in the ABL and hence the available wind resources for energy production. Figures 4(b) and 4(f) show the magnitudes of the planarand temporally averaged vertical kinematic momentum flux τ = u w 2 + v w 2 with u w = uw + τ xz − u · w . In the barotropic cases, the vertical kinematic momentum flux has its maximum at the ground and monotonically decreases with increasing height, approaching zero aloft where the influence of the surface is negligible. In contrast, shear is not only generated at the surface, but also due to the changes in geostrophic wind in the baroclinic cases. Right above the LLJ, where dG/dz < 0, a local maximum of momentum flux forms. Such local maxima have been reported previously by Conangla and Cuxart [20], who concluded that the geostrophic shear increases the wind shear and turbulence production away from the surface. The flow above the LLJ is thus shear dominated. In the neutral cases, the differences between the baroclinic and barotropic cases are less pronounced than in the SBL cases. This could result from the enhanced turbulence (see Table I) and the lower imposed baroclinicity strength in the neutral cases.

III. BOUNDARY-LAYER CHARACTERISTICS
To estimate the baroclinicity strength S 0 , we follow the definition of Momen et al. [24]. For negative shear baroclinicity, where the geostrophic wind G does not change its direction with height, the baroclinicity strength is defined as Momen et al. [24] defined the baroclinicity strength based on the gradient measured at the surface. However, as dG/dz = 0 at the surface in the SBL cases, we define S 0 based on (dG/dz) max . Following Momen et al. [24], we approximate the ABL height z i by the base height of the capping inversion layer, which is z i = 1200 m for all cases.
The values for S 0 are listed in Table II. Thus, the stable cases have a stronger baroclinicity than the neutral cases, while the baroclinicity is strongest in the baroclinic low cases. The reason is that the gradient of the geostrophic wind compared to its absolute values is highest for the baroclinic low cases. As a result, the local maximum in the total vertical momentum flux, which arises due to baroclinicity, is most pronounced for the SBL baroclinic low case, where the baroclinicity strength is highest [see Figs. 4(b) and 4(f)].
The planar averaged turbulence intensity at hub height TI = 1/3(u 2 + v 2 + w 2 ) / u 2 + v 2 + w 2 , with u 2 = u 2 + τ xx − u 2 , is given in Table II. The baroclinic low cases show a reduced or equal turbulence intensity in comparison to their respective barotropic cases, for which G 0 is the same. An increase in TI is observed for the baroclinic high cases. In agreement with the time-and horizontalaveraged hub height velocity at the inlet V inflow , the friction velocity and the Obukhov length are lowest for the baroclinic low cases and highest for the baroclinic high cases (see Table II). The planar-and temporally averaged potential temperature θ is shown in Figs. 4(c) and 4(g). For the SBL, the surface is cooled such that a stable stratification exists near the surface, followed by a neutral temperature profile up to the capping inversion. Above the capping inversion, there is stable free-atmospheric stratification. In the NBL, the potential temperature is constant below the capping inversion, above which the atmosphere is stably stratified. There are only minor differences in the potential temperature profiles in barotropic and baroclinic cases.
Because of the imposed geostrophic shear, the Ekman spiral is altered, yielding different wind angles for different cases. Figures 4(d) and 4(h) show the mean wind angle α = tan −1 v / u . Note that, due to the applied windangle controller, the wind direction has a constant value of α = 0 • at hub height for all cases under consideration. Typically, the wind veer is more pronounced in a SBL than in a NBL, and stronger for cases with negative geostrophic shear. A changing wind veer can alter the spatial structure of the wakes [56] and is, consequently, important to take into consideration.

IV. EFFECT OF BAROCLINICITY ON FLOW IN AND AROUND WIND FARMS
In the previous section, we showed that negative geostrophic shear alters the mean wind profile and causes a lower jet height. Furthermore, above the LLJ, additional shear and turbulence are created by baroclinicity. In the following, we examine how these changes affect the flow in and around a wind farm. Figure 5 shows the instantaneous velocity at hub height (x-y plane) for the (a) SBL barotropic, (b) SBL baroclinic low, and (c) NBL baroclinic low cases. The wakes behind the first row for the SBL barotropic and SBL baroclinic low cases are straight and have a similar velocity deficit. In contrast, the wakes in the second row seem to be shorter and meander stronger in the SBL baroclinic low than in the SBL barotropic case. Further downstream, the velocity TABLE II. The columns from left to right indicate the case name, the LLJ height (z LLJ ), the velocity at the LLJ height (v LLJ ), velocity (V inflow ) and turbulence intensity (TI) at hub height, frictional velocity (u * ), Obukhov length (L), and baroclinicity strength (S 0 ).  deficit in the SBL baroclinic low case is higher and persists longer than in the SBL barotropic case. This difference is caused by the height of the LLJ, which is located above the turbines in the SBL barotropic case, while it overlaps with the rotor swept area in the SBL baroclinic low case. In the latter case, the first rows extract most of the energy from the LLJ, leaving fewer resources for downstream turbines. Furthermore, we observe a significant clockwise deflection of the wakes downstream. This is likely a result of downstream entrainment of momentum from the layer above [79,80], where the flow direction is rotated [see Fig. 4(d)].
In the NBL case [ Fig. 5(c)], the atmospheric scales are larger and form high-and low-velocity streaks. Additionally, the turbulence intensity is higher in the NBL cases; see Table II. Consequently, the meandering of the wakes in the lateral direction is more pronounced. Figures 5(d)-5(f) show y-z planes of the instantaneous velocity at a distance of 1D behind the last turbine row. In the SBL barotropic case [ Fig. 5(d)] the flow above the LLJ is nonturbulent. However, in both baroclinic cases [Figs. 5(e) and 5(f)] the flow is more turbulent above the jet due to the geostrophic shear. Furthermore, due to the prevailing wind veer, the wakes are skewed significantly in the SBL cases [56]. In the NBL, this effect is less pronounced as the veer is much smaller. The jet in the NBL case is stretched up to a height of 1.2 km, and this layer is well mixed and more turbulent than in the SBL cases. For a further discussion of the effect of baroclinicity on the flow structures, we refer the reader to Momen et al. [24]. Figure 6 shows the time-averaged velocity, averaged over the spanwise extent of the wind farm, which clearly reveals the LLJ height. The LLJ is lowest in the SBL baroclinic low case [ Fig. 6(b)] in which it resides largely within the rotor swept area of the turbines. When the LLJ occurs roughly at hub height, the first rows of the wind farm extract most of the momentum from the LLJ. Consequently, the energy that is available for downstream turbines is significantly reduced [36]. Figure 6(b) shows that the jet velocity is reduced towards the end of the farm. The effect is much less pronounced in the other cases, where the LLJ is positioned higher.
To show how this phenomenon affects the velocity at hub height, we evaluate the velocity averaged over the spanwise extent of the wind farm, normalized by the inflow velocity; see Fig. 6(d). For all cases, the velocity reduces significantly at the turbine locations, and there is a decreasing trend in the downstream direction over the length of the wind farm. Compared to all other cases, the SBL baroclinic low case exhibits a greatly reduced downstream wake recovery. Furthermore, we observe that the wake behind the first turbine row recovers faster for the baroclinic cases, while the barotropic cases show increased wake recovery further downstream in the wind farm. Finally, the wake recovery is significantly slower in stable than in neutral conditions, in agreement with previous observations by Abkar and Porté-Agel [48]. The wake recovery further downstream in the wind farm depends strongly on the entrainment of energy from the flow above the wind farm [1]. Therefore, we analyze the entrainment fluxes to understand the slower wake recovery for baroclinic cases in the following section.

V. ENERGY BUDGET ANALYSIS
To study how negative geostrophic wind shear affects entrainment fluxes and the power production of the turbines in the wind farm, we perform an energy budget analysis [81].
To obtain the total kinetic energy equation, the momentum equation (2) is first multiplied with u i and then averaged over time. Subsequently, we integrate the resulting equation over a control volume to isolate the energy balance around each turbine row. We select control volumes ∀ of size s x × 5s y × D surrounding the center of each turbine row. The corresponding energy equation is given by Here, P is the power production by a turbine row; E k is the mean-flow transport of kinetic energy, including the meanflow and resolved turbulent kinetic energy; T t and T SGS are the transport of momentum by resolved and SGS turbulent fluxes, respectively. The transfer of energy due to pressure fluctuations is represented by F, while B is the turbulence destruction or production due to buoyancy, and G is the mean geostrophic forcing. Figure 7 shows the energy budget analysis for the (a) SBL barotropic and (b) SBL baroclinic low cases, normalized by the power produced by the first row of turbines. 023007-9  Table I. All the terms are normalized by the power production of the first turbine row. The symbols in the legend are defined in Eq. (9).
For brevity, this analysis is only shown for the two cases that reveal the effect most pronounceable. The trends for the other cases are similar to the SBL barotropic case. We denote by P, D, G, and B the energy sinks, and by E k and T t the energy sources. Using an energy budget analysis for wind farms with an infinite lateral size in barotropic NBLs, Meyers [81] and Cortina et al. [44] have shown that the kinetic energy contribution E k is highest at the first row and shows a decreasing trend over the length of the farm since the turbines extract kinetic energy from the flow. Meanwhile, the turbulent transport T t increases as turbulent wakes form and interact behind the turbines. In Fig. 7 we observe the same trend and show that, for the SBL barotropic case, the turbulent transport T t becomes the dominant contribution, rather than E k , starting from the fourth row. However, this transition takes place further downstream for the SBL baroclinic low case because the lower jet height limits the energy entrainment into the wind farm. Abkar and Porté-Agel [82] and Cortina et al. [44] have shown that the transition point, at which T t starts exceeding E k , depends on the wind farm layout, while Abkar and Porté-Agel [82], Cortina et al. [83], Wu and Porté-Agel [50], and Gadde and Stevens [64] further showed that the transition point also depends on the atmospheric stratification. In particular, T t is higher for denser wind farms [44,82,84] and for lower atmospheric stability [50,64,82,83]. Furthermore, Gadde and Stevens [36] showed that T t is highly dependent on the height of the low-level jet. In the SBL baroclinic low case, the lower jet height limits the energy entrainment into the wind farm and consequently the vertical energy transport from above only becomes dominant further downstream in the wind farm.
To demonstrate the reduced energy entrainment from above for the SBL baroclinic low case, we show the integrated vertical entrainment flux through the top (z h + D/2) and bottom (z h − D/2) planes of the control volume in Fig. 8. In both SBL and NBL cases, the barotropic cases exhibit more entrainment from above than the baroclinic cases. The entrainment is smallest for the SBL baroclinic low case, where the LLJ is lowest and resides within the rotor swept area. The entrainment flux from below is much smaller than from above and almost unaffected by baroclinicity.
The decreased downward flux in the SBL baroclinic cases is also visualized in Fig. 9. The reason for the lower entrainment from above in the baroclinic cases is twofold. (i) The velocity above the turbines, relative to the velocity at hub height, is higher for the barotropic case than for the baroclinic cases (see Fig. 10 in which we normalize the velocity by the hubheight velocity, the representative for power production of the first turbine row). Consequently, there is relatively less momentum available for downward entrainment in baroclinic cases. (ii) The negative shear in the baroclinic cases gives rise to upward turbulent flux, which reduces the net energy entrainment into the wind farm. The reason is that the turbulence and shear profiles above the nose of the LLJ are different in the barotropic and baroclinic conditions; see Figs. 4(a), 4(b), 4(e) and 4(f). The baroclinic cases have a larger negative shear above the wind maximum than the barotropic cases. The term u w ∂u/∂z in the energy budget has to be negative to produce turbulence. For the baroclinic cases, there is a large negative shear above the nose of the jet, i.e., ∂u/∂z < 0, which results in a positive vertical entrainment flux aloft, i.e., u w > 0.
This upwards (positive) entrainment flux due to the negative shear is clearly visible in Figs. 9(d)-9(f). The highest upwards flux outside the IBL of the wind farm is observed for the SBL baroclinic low case, while the vertical entrainment outside the IBL of the wind farm is negligible in the SBL barotropic case.
In summary, the presence of a wind farm in the ABL creates a momentum deficit region and therefore causes a downward entrainment flux in the presence of positive mean shear in the ABL [38,[84][85][86]. However, the large negative shear created by the baroclinicity in the SBL counteracts this effect and therefore reduces the net entrainment downwards. The upward turbulent flux changes the momentum transfer from the jet and hinders downward energy entrainment into the wind farm. We note that none of the cases considered in this study reaches a fully developed wind farm regime, where the statistical properties of the flow are constant along the streamwise direction such that the energy available at the wind turbine locations must be entrained entirely from the layers above [38,84,85,[87][88][89]. In contrast to infinitely wide wind farms considered in previous studies [38,52,84,[87][88][89], wind farms with a finite spanwise size (see also, e.g., Refs. [41,47,64,90,91]) are expected to be influenced by entrainment from the sides due to the prevailing wind veer. Wu and Porté-Agel [50] and Allaerts and Meyers [81] assumed that a fully developed region might not be reached in various atmospheric conditions even for large wind farms. However, Cortina et al. [44] stated that the flow physics are very similar to the physics in fully developed regimes, once T t dominates. Figure 11 shows the power production normalized by the power production of its first row for the barotropic reference case. We observe that the total power production strongly depends on the baroclinicity. To be specific, the baroclinic high cases have the same geostrophic wind at high elevations as the barotropic cases. Here, baroclinicity leads to higher wind velocities at hub height and consequently the highest power production. The SBL barotropic case and the SBL baroclinic low case have similar velocities at hub height, as shown in Table II. However, the trend in power production is quite different, which will be discussed below. Figure 11 highlights that baroclinicity alters the mean wind profile in the ABL and thereby affects the wind energy available for extraction by the turbines. Consequently, including baroclinicity effects in numerical and analytical models that are used to design wind farms is crucial to obtaining accurate performance predictions. Figure 12 shows the wind farm power production, normalized by the power production of its first turbine row for the (a) SBL and (b) NBL cases. Behind the first row, the wind turbine wakes cause a strong drop in power production. The SBL baroclinic low case has the highest normalized power production for the second and third rows. This is caused by the high kinetic energy of the LLJ, which resides at turbine rotor height in this case [ Fig. 7(b)]. Further downstream, the power production of the SBL baroclinic low case is lowest due to the previously discussed lack of energy entrainment from above. The normalized power production of the SBL baroclinic high case is slightly below the reference barotropic case. For the NBL cases [ Fig. 12(b)], the effect of the negative geostrophic wind shear on the normalized power production is negligible. This is because the wind profile in the region directly above the wind farm is similar in the barotropic and baroclinic cases; see Fig. 10.

Turbine row
Turbine row Relevant differences in the available momentum above the wind farm are only observed above approximately 400 m. Consequently, the entrainment into the wind farm is only slightly lower for the baroclinic cases than for the barotropic case (Fig. 8).

VII. CONCLUSION
We perform LESs to study the effects of negative shear baroclinicity on wind farm performance in stable and neutral ABLs. We find that even modest negative geostrophic wind shear of 3 ms −1 km −1 with height can greatly alter wind farm power production, as geostrophic wind shear modifies the mean wind profile in the ABL significantly [see Figs. 3(a), 3(e), and 10]. For the baroclinic cases considered in this study, which have the same geostrophic wind velocity at either the surface or the top of the domain as the barotropic cases, the velocities at hub height differ by up to 2.3 m/s. An increase in the geostrophic wind near the surface relative to barotropic conditions can significantly increase the velocity at hub height and, consequently, the absolute power production of a wind farm (the baroclinic high cases). Meanwhile, a reduction of the geostrophic wind aloft with respect to barotropic conditions can lower the LLJ to turbine height (the SBL baroclinic low case). In that case, normalized power production increases for the first turbine rows, but decreases significantly for downstream rows, as less momentum is available above the wind farm for downward entrainment.
The negative geostrophic wind shear does not only alter the mean velocity profile, but also modifies turbulence in the ABL. To be specific, the negative geostrophic wind shear increases turbulence and vertical momentum fluxes at higher altitudes. Typically, additional turbulence aloft is considered beneficial for the performance of extended wind farms as it can aid wake recovery and entrainment. However, an energy budget analysis reveals that negative shear baroclinicity reduces entrainment into the wind farm. In fact, the negative shear creates a positive turbulent flux above the LLJ and causes an upward transport of momentum. This limits the entrainment of energy from the jet into the wind farm and ultimately hinders the performance of turbines further downstream in the wind farm.
The observed alterations of wind resources and turbulence in the ABL and resulting changes in wind farm power production indicate that baroclinicity should be considered in wind farm modeling and power production forecasts. Therefore, future studies are required to incorporate the effects of negative shear baroclinicity in analytical models. Furthermore, it will be important to investigate different forms of baroclinicity [24,27], resulting from different alignment between pressure and temperature gradients. This includes positive shear, as well as thermal advection. Besides, it is important to note that there are no physical constraints on the relative orientation between T and P, and T might vary with height [24]. Consequently, it will be relevant to investigate the impact of nonlinear changes in the magnitude of the geostrophic wind velocity with height and to consider time-dependent geostrophic wind.

ACKNOWLEDGMENTS
This work is part of the Shell-NWO/FOM initiative Computational Sciences for Energy Research of Shell and Chemical Sciences, Earth and Live Sciences, Physical Sciences, FOM, and STW. R.J.A.M.S. gratefully acknowledges funding by STW VIDI, under Grant No. 14868. This work was partly carried out on the Dutch national einfrastructure with the support of SURF Cooperative. We acknowledge PRACE for awarding us access to MareNostrum at the Barcelona Supercomputing Center (BSC), Spain (projects 2020225335 and 2020235589).

APPENDIX
Here, we show the time-averaged version of Fig. 5 to give a reference flow field in Fig. 13. The figure shows the time-averaged velocity at hub height for the (a) SBL barotropic, (b) SBL baroclinic low, and (c) NBL baroclinic low cases. Starting from the third turbine row, the velocity deficit in the SBL baroclinic low case is higher than in the SBL barotropic case. In addition to the phenomena discussed for the instantaneous flow field in Fig. 5, the turning of the wakes is clearly visible in both stable cases. The strongest wake turning (5 • ) is observed in the last turbine row for the SBL baroclinic low cases. For the NBL baroclinic low case, the wind veer is less, namely 1.5 • in front of the last row.
To verify that the fixed orientation of the rotor heads perpendicular to the x direction in our simulations do not affect our findings, we perform simulations on a coarse grid with a resolution of x = 24 m, y = 15 m, z = 5 m using actuator disks, for the stable cases. Figure 14 shows the power production for the stable cases: (1) the rotor heads are oriented perpendicular to the x direction and (2) the rotor heads are oriented perpendicular to the local flow angle. The difference in total power production between the two cases is small ( P tot, local angle / P tot, x aligned ≤ 2.5%) and, most importantly, the same trends are observed when comparing baroclinic and barotropic cases.