Interaction between low-level jets and wind farms in a stable atmospheric boundary layer

Low-level jets (LLJs) are the wind maxima in a weak or moderately stable atmospheric boundary layer. The momentum in the LLJs can serve as a source of high capacity energy for wind turbines. In this work, large-eddy simulations are used to study the effect of LLJs and low-lying temperature inversions on the power production of wind farms by systematically increasing the stability of the boundary layer. The LLJs increase the shear in the boundary layer and enhance power production. Flow visualization reveals that the low-lying capping inversion limits the growth of the wind farm internal boundary layer and forces the wind to go around the wind farm. This effect is stronger at higher stratification; i.e. more wind goes around the wind farm. An energy budget analysis is carried out to analyze the effect of different flow phenomena on turbine power production. The budget analysis shows elevated entrainment rates in the presence of LLJs, emphasizing the effect of jets as an additional momentum source. Furthermore, the height of the LLJ determines the wake recovery and power production in the wind farm. When the jets are far above the wind turbines, higher entrainment in the rear of the wind farm causes an increase in power production. If the jets are at the same height as the top of the wind turbines, the wind turbine wake meandering in the vertical direction can directly extract the energy in the jets. At high stratification, the combined effect of buoyancy destruction and turbulence dissipation is more than the turbulent entrainment, and causes a slower wake recovery.


Introduction
The atmospheric boundary layer is dynamic and undergoes continuous transitions during the day due to the changes in forcing, such as the surface heat flux and the geostrophic wind. During the evening boundary-layer transition with cooling at the ground, the flow above the surface layer decouples from the surface friction. Consequently, the balance between the Coriolis, frictional, and pressure forces is disturbed, and the flow above the surface layer accelerates. The acceleration produces a super-geostrophic jet at the top of the nocturnal stable boundary layer (SBL) at low heights (50 -1000 m) of the atmosphere (Smedman et al. 1996). This super-geostrophic wind is called the low-level jet (LLJ), and it generally forms due to the aforementioned frictional decoupling combined with inertial oscillations (Blackadar 1957;Thorpe & Guymer 1977). Additionally, a LLJ can also form due to large-scale baroclinicity or the pressure gradient due to cooling over sloped terrains (Mahrt 1999). In general, LLJs form in nocturnal conditions when there is surface cooling. Mahrt (1998) classified the nocturnal boundary layer into three stability regimes: 1) the weakly stable regime, characterized by continuous turbulence and a small downward heat flux, which is limited by the temperature fluctuations, 2) transition stability regime, where the quantities change rapidly with the increasing stability and the downward heat flux reaches a maximum, and 3) the very stable regime where the downward heat flux is small, limited by the turbulent vertical fluctuations, which are suppressed by buoyancy. LLJs of practical importance are characterized by high shear and weak to moderate stability (Baas et al. 2009;Banta 2008). The shear in the LLJ is strong enough to generate continuous turbulence with maximum and minimum values at the surface and top of the SBL, respectively (Mahrt 1998).
LLJs are frequently observed in many parts of the world, with occurrences in the Western ghats of India (Prabha et al. 2011), the Great Plains of the United States (Kelley et al. 2004;Banta et al. 2002;Lundquist 2003), the North sea, and the Baltic sea of Europe (Kalverla et al. 2019;Smedman et al. 1993). LLJs have wide-ranging impact in the atmosphere, the wind shear in an LLJ can cause dangerous deviations from the nominal glide path of an aircraft resulting in a premature touch-down (Wittich et al. 1986), the high-speed jet can increase the dispersion and transport of pollutants in the atmosphere (Beyrich 1994), and a LLJ formed in the upper layers of the atmosphere can even impact bird migrations (Liechti & Schaller 1999). Furthermore, LLJs have a direct impact on wind energy-based power production (Sisterson & Frenzen 1978). It is a common practice in wind power assessment to use simple, power-law velocity profiles. However, this overlooks the high power density and strong shears associated with a LLJ phenomenon. This results in an underestimation of both power production and the fatigue loads on the wind turbine (Gutierrez et al. 2017). Wilczak et al. (2015) report that the LLJs increase the capacity factors by over 60% under nocturnal conditions. Present-day wind turbines are reaching heights above 200 m due to which interactions with LLJs becomes unavoidable. Consequently, it is imperative to study the interaction between LLJs and wind farms.
When a large number of wind turbines operate in a wind farm, the structure of the boundary layer changes due to the momentum extraction by the turbines. Both numerical simulations and wind tunnel experiments show the development of an internal boundary layer (IBL) at the entrance of the wind farm (Frandsen et al. 2006;Calaf et al. 2010;Chamorro & Porté-Agel 2011). Further downstream, in the fully developed regime, all the momentum is derived from vertical entrainment (Calaf et al. 2010;Cal et al. 2010). Owing to the simplicity, most wind farm and atmospheric boundary layer simulations in the past have focused on pressure-driven neutral boundary layers. The underlying assumption in such simulations is that the wind turbines reside in the inner regions of the atmospheric boundary layer, where the outer layer effects such as the rotation of the Earth and thermal stratification are negligible (Stevens & Meneveau 2017). However, the wake recovery and entrainment of fresh momentum from outside the IBL are heavily dependent on the turbulence intensity, which in turn depends on the stratification in the atmosphere (Abkar & Porté-Agel 2015; Barthelmie et al. 2015). Furthermore, the wind follows an Ekman spiral due to the Coriolis force and causes a wind veer affecting the turbine wakes as well as the combined wind farm wake. In essence, neglecting the stratification and Coriolis forces is too simplistic when large wind farms are considered.
In the past decade, the large-eddy simulation (LES) technique has been extensively used to study turbulence in the atmosphere (Moeng 1984;Mason & Thomson 1992), and the interaction between the atmospheric boundary layer and wind farms (Stevens & Meneveau 2017;Stevens et al. 2014b;. LES has been successfully used to simulate both convective and stable atmospheric boundary layers at both weak and moderate stratification (Mason 1989;Nieuwstadt et al. 1993;Mason & Derbyshire 1990;Saiki et al. 2000;Kosović & Curry 2000). Numerical simulations of weak and moderately stratified SBL are easier because of the continuous turbulence and the absence of global intermittency (Mahrt 2014). The simulations of highly stratified boundary layers are challenging due to the mesoscale motions, gravity waves, the unsteady nature of the boundary layers, and LLJs. Nocturnal LLJs under weak to moderate stratification have been extensively studied with LES (Beare et al. 2006;Kosović & Curry 2000).
Recently, Allaerts & Meyers (2017) studied the impact of the capping inversion on the power production of 'infinitely' wide wind farms in conventionally neutral boundary layers. They report that the IBL pushes the capping inversion upwards, which generates pressure perturbations that travel upstream as gravity waves and slow down the upstream wind. Furthermore, Allaerts & Meyers (2018) studied wind farms in a neutral-to-stable boundary layer transition, and they note that in a steady-state SBL, the LLJ impact the power production. Also, measurements and LES studies of wind farms in SBL (Allaerts & Meyers 2018;Dörenkämper et al. 2015;Witha et al. 2014;Barthelmie et al. 2015) show that due to low turbulence intensity, wake recovery is reduced compared to the unstable and neutrally stratified boundary layer. In addition, the rotation of the Earth affects the power production with the Coriolis force rotating the wind farm wake (van der Laan & Sørensen 2017). For certain wind directions it has been found that even the horizontal component of the Earth's rotation influences the turbulent fluxes in a wind farm (Howland et al. 2020). Furthermore, the vertical wind veer in the Ekman spiral causes a skewed spatial structure of the turbine wake and also enhances shear production of turbulent kinetic energy leading to larger flow entrainment and thus faster wake recovery . It is a common practice in the wind energy community to use periodic boundary conditions in the spanwise direction, which results in 'infinitely' wide wind farms (Calaf et al. 2010;Stevens et al. 2014b;Allaerts & Meyers 2015;Wu & Porté-Agel 2013). However, in the presence of Coriolis force, which induces appreciable spanwise flow, this assumption might lead to under-prediction of power in the outer turbine rows.
Previously, wind turbine and LLJ interactions have been studied by Lu & Porté-Agel (2011), who performed LES of the flow over a turbine in a doubly periodic domain (an 'infinite' wind turbine array) with actuator line modeling, and they report non-  Figure 1: Sketch of the essential flow phenomena in wind farms in a SBL, including wakes and their superposition, the entrainment of momentum from above, and the development of the IBL. On the left, the typical temperature and velocity profiles, which reveal the LLJ and the capping inversion, are sketched. axisymmetric turbine wakes and LLJ elimination due to energy extraction by the turbines. Furthermore, with LES of an infinite wind farm in a diurnal cycle and observed the formation of LLJ in nocturnal conditions ; Shamsoddin & Porté-Agel (2017) also report that the LLJ is weakened because of the extraction of momentum by the turbines. Similar results have been reported by Fitch et al. (2013) in the mesoscale weather model simulations of an infinite wind farm. Recently, Na et al. (2018) performed LES of a small wind farm with 12 turbines arranged in 3 rows and 4 columns with a LLJ above it in a spanwise periodic domain. They report faster wake recovery due to the enhanced vertical kinetic energy flux created by the LLJ. The studies mentioned above do not provide a complete picture of the interaction between LLJs and wind farms, either because of the 'infinite' wind farm assumption, or the spanwise periodic small size wind farms. Furthermore, these studies of LLJ and wind farm interaction have not focused on the flow adjustment in and around the wind farm and the impact of different jet heights on wind farm power production. The flow is further complicated by interaction between the low-lying capping inversion and the wind farm IBL. It is necessary to understand the coupling between stable stratification, flowadjustment, and LLJs on wind farm power production. Figure 1 shows the essential flow phenomena such as the IBL growth, turbine wake recovery, modification of the capping inversion, and the entrainment of momentum from above by turbulence.
In this work, we study the effect of LLJs in stable stratification in a finite wind farm (with no periodicity along the boundaries) by systematically varying the surface cooling rate. To the best of our knowledge, this has not been studied earlier. The objective of the study is two-fold, first, to study the effect of LLJ and stable stratification on the power production of a wind farm, and second, to study the effect of low-lying capping inversions on the flow adjustment in and around a 'finite' wind farm.
The remainder of the paper is structured as follows, in § 2, the LES method and numerical method are detailed. In § 3 important boundary layer properties, the IBL growth above the wind farm, and the flow adjustment around the wind farm are discussed. In § 4, an analysis of the different flow phenomena is carried out by performing an energy budget analysis. Furthermore, in § 5, the finite wind farm effect due to the wind veer is discussed, which is followed by the conclusions in § 6.

Large-eddy simulations
In LES, the flow features larger than the filter size are fully resolved, while the sub-filter size eddies are modeled. Our code is an updated version of the one used by (Albertson & Parlange 1999;Calaf et al. 2010). This updated code has been validated for neutral and SBLs as well as the flow-through wind farms Zhang et al. 2019;Gadde & Stevens 2019). The governing equations and numerical method are discussed in § 2.1, and the boundary layer initialization is explained in § 2.2, followed by the wind farm setup in § 2.3.

Governing equations and numerical method
The LES code we use integrates the filtered Navier-Stokes equations written for a wallbounded turbulent flow (Albertson 1996) and employs the Boussinesq approximation to model buoyancy. The f−plane approximation is used to model the Coriolis forces. The governing equations are: where, the tilde represents spatial filtering, u i = ( u, v, w) and θ are the filtered velocity and potential temperature, respectively, g is the acceleration due to gravity, β = 1/θ 0 is the buoyancy parameter with respect to the reference potential temperature θ 0 , δ ij is the Kronecker delta, f c is the Coriolis parameter. The boundary layer is driven by a mean pressure gradient, p ∞ , represented by the geostrophic wind with, U g = − 1 ρfc ∂p∞ ∂y and V g = 1 ρfc ∂p∞ ∂x as its components. p = p * /ρ + σ kk /3, is the modified pressure obtained by adding the trace of the sub-filter scale stress, σ kk /3, to the kinematic pressure or pressure perturbation, p * /ρ, where ρ is the density of the fluid. f i = ( f x , f y , 0) represents the turbine forces, which are modeled using a filtered actuator disk approach (á. Jiménez et al. 2010;Calaf et al. 2010Calaf et al. , 2011. The molecular viscosity is neglected as it is a high Reynolds number flow, which is a common practice in atmospheric boundary layer simulations. τ ij = u i u j − u i u j is the traceless part of the sub-filter scale stress tensor, and q j = u j θ − u j θ is the sub-filter scale heat flux tensor. The sub-filter stresses and heat fluxes are modeled as, where, S ij = 1 2 (∂ j u i + ∂ i u j ) is the filtered strain rate tensor, ν T is the eddy viscosity, C s is the Smagorinsky coefficient for the sub-filter stresses, ∆ is the filter size, ν θ is the eddy heat diffusivity, D s is the Smagorinsky coefficient for the sub-filter scale heat flux, and | S| = 2 S ij S ij . We use a tuning-free, scale-dependent model based on Lagrangian averaging of the coefficients (? Stoll & Porté-Agel 2006, 2008 to dynamically calculate the Smagorinsky coefficient. The error in the calculation of the Smagorinsky coefficients is minimized over fluid pathlines preserving the local fluctuations of the coefficients. This makes the model particularly suitable for inhomogeneous flows, such as the flow through a wind farm or over complex terrain. We use a pseudo-spectral method to calculate the partial derivatives in the streamwise and spanwise directions. The vertical direction is treated with a second-order central difference method. The solution is advanced in time by a second-order accurate Adams-Bashforth scheme. The aliasing errors resulting from the folding back of high wavenumber energy to the resolved scales due to the calculation of non-linear terms in physical space is prevented by using a 3/2 anti-aliasing method (Canuto et al. 1998). For pointwise energy conservation, the convective term in the equation 2.2 is written in the rotational form (Ferziger & Perić 2002). More information about the numerical method can be found in Albertson (1996). The computational domain is discretized uniformly with n x , n y , and n z points in the streamwise, spanwise, and vertical directions, respectively. Therefore, the corresponding grid sizes are ∆ x = L x /n x , ∆ y = L y /n y , and ∆ z = L z /n z , where L x , L y , and L z are the dimensions of the computational domain. The computational grid is staggered in the vertical direction with the first grid point for u, v, and θ located at a distance ∆ z /2 above the ground. The computational plane for the vertical velocity, w, is located at the ground. No-slip and free-slip boundary conditions with zero vertical velocity, w = 0, are used at the top and bottom boundaries, respectively. In wall-modeled LES of atmospheric boundary layers, the first grid point generally lies in the surface layer and the Monin-Obukhov similarity theory (Moeng 1984) can be used to model the instantaneous shear stress τ i3|w and buoyancy flux q * at the wall as follows and where u i and θ represents the filtered velocities and potential temperature at the first grid point respectively, u * is the frictional velocity, z o is the roughness length, κ is the von Kármán constant, u r = √ u 2 + v 2 is filtered velocity magnitude at the first grid level, and θ s is the filtered potential temperature at the surface. ψ M and ψ H are the stability corrections for momentum and heat flux, respectively. For SBLs, we use the stability correction suggested by Beare et al. (2006), i.e. ψ M = −4.8z/L and ψ H = −7.8z/L, where L = −(u * 3 θ 0 )/(κgq * ) is the surface Obukhov length.

Boundary layer initialization
Simulating strongly stratified boundary layers with LES is complicated due to the presence of globally intermittent turbulence (Mahrt 2014). In the present work, we consider a moderate stratification (z i /L ≈ 2, where z i is the boundary layer height), which remains continuously turbulent. The boundary layer represents a typical quasiequilibrium moderately SBL, with a pronounced LLJ similar to those observed over polar regions and equilibrium night-time conditions over land in mid-latitudes. The case is well-documented under the global energy and water cycle experiment atmospheric boundary layer study (GABLS−1) initiative and LES inter-comparison studies (Kosović & Curry 2000;Beare et al. 2006). The initial potential temperature profile has a mixed layer (with constant potential temperature 265 K) up to 100 m with an overlying capping inversion of strength 0.01 Km −1 . The reference potential temperature and roughness length are set to 263.5 K and 0.1 m, respectively, and a constant surface cooling is applied. The boundary layer is driven by the geostrophic wind with the horizontal components G = (U g , V g ) = (8.0, 0.0) ms −1 . The Coriolis parameter is set to f c = 1.39 × 10 −4 s −1 (corresponding to latitude 73 • N). The initial wind profile is set equal to the geostrophic wind except. To spin up turbulence, uniformly distributed random perturbations with an amplitude 3% of the geostrophic wind are added to velocities below a height of 50 m. Similarly, uniformly distributed random perturbations of magnitude 0.1 K are added to the initial temperature profile. Detailed information about the SBL can be found in Beare et al. (2006). The boundary layer reaches a quasi-steady state at the end of 8 th hour. The quasi-steadystate is said to have been reached when the temperature profile changes at a constant rate while the velocity and other turbulent quantities have reached a steady-state (Kosović & Curry 2000). Our code has been validated for the GABLS-1 boundary layer with a cooling rate of 0.25 K · hour −1 . Extensive grid resolution studies were carried out, and it was found that the Lagrangian averaged scale-dependent (LASD) model produces accurate results when compared with the Smagorinsky model, which agrees with the previous study by Stoll & Porté-Agel (2008). We performed five simulations with surface cooling rates C r = [0.0, 0.125, 0.25, 0.375, 0.5] K · hour −1 . Details about the different cases are documented in table 1.

Wind farm setup
We consider a large wind farm with 40 wind turbines. The turbines are distributed in an array of 4 rows and 10 columns. The turbine diameter is D = 90 m and the hub height is z h = 90 m. The turbines are separated by a distance of s x = 7D and s y = 5D in the streamwise and spanwise directions, respectively. The computational domain is 11.52 km × 4.6 km × 3.84 km. The details of the computational domain and wind farm layout are given in figure 2. We use the Monin-Obukhov similarity theory to model the wall stresses, according to which the first grid point above the ground should be in the inertial sublayer. For the GABLS-1 case, cautioning against using very high resolution near the ground which violates the similarity theory, Basu & Lacser (2017) prescribe a guideline of z 1 50z o , where z 1 the height of the first grid point above the ground. Accordingly, we fix the vertical grid resolution to be 5 m. We performed extensive grid resolution studies with this vertical resolution and found that a resolution of 9 m resolves all the relevant horizontal scales. This amounts to a resolution of 9 m × 9 m × 5 m which is found to be sufficient to resolve all the relevant length scales. The domain is discretized by 1280 × 512 × 768 grid points in the streamwise, spanwise, and vertical directions, respectively. The computational domain has approximately 500 million grid points. We use a large vertical extent and a Rayleigh damping layer (Klemp & Lilly 1978) with a strength of 0.016 s −1 in the top 25% of the domain to reduce the effects of gravity waves. It was found that this damps out most gravity waves that are generated. To ensure that the streamwise fringe layer does not affect the turbulence statistics, we performed a simulation in a domain of size 17.28 km × 4.6 km × 3.84 km. The streamwise fringe layer in this domain is at a distance 75D away from the wind farm. We find that the streamwise domain size does not affect the turbulence statistics relevant to the study.
To obtain realistic inflow conditions, we employ the concurrent precursor technique (Stevens et al. 2014a). In this technique, simulations are run in two domains concurrently. Firstly in the precursor domain, we perform the atmospheric boundary layer simulations without wind turbines to generate inflow conditions, and secondly, the quantities from the precursor domain are used as the inlet conditions for the wind farm domain. The forcing is done in the wind farm domain by gradually blending the velocities in the fringe layer. Due to the Coriolis forces, an Ekman spiral, which induces considerable spanwise flow, is formed. Therefore, we use fringe layers in both the streamwise and spanwise direction to eliminate the effects of the periodic boundaries. We fix the fringe layer length to be 10% of the computational domain in the streamwise and spanwise directions, respectively. The equilibrium wind angle under geostrophic forcing depends on the stability conditions, which results in a different geometric pattern of the turbines and complicates the analysis. To ensure the same farm layout in all simulations, we use a proportional-integral (PI) controller (Sescu & Meneveau 2014), similar to the one used by Allaerts & Meyers (2015), to rotate the incoming flow such that the planar averaged wind angle at hub height is always zero. Even then, local changes in the wind angle upstream of a turbine can result in turbine yaw misalignment, which results in sub-optimal energy production. To prevent this, each turbine in the simulations has an individual yaw-angle controller, which reorients the orientation of the turbines such that it is perpendicular to the local wind angle measured 1D upstream of each turbine.

LES of a finite wind farm
All the simulations are carried out in two stages. In the initial or the spin-up stage, only the SBL in the precursor domain. The SBL reaches quasi-steady state, at the end of the 8 th hour. In the second stage, the turbines are introduced in the main domain, and the simulation in both domains are continued concurrently for one more hour in which the transient effects of the turbine startup subside. Finally, both the simulations are run for one more hour, and the statistics are collected in the last hour i.e. the 10 th hour. The one-hour statistics was collected in six ten minute samples to quantify the uncertainty in the mean values. The standard deviation over the mean is found to be within 5% for all the equilibrium profiles. In the following sections, the basic boundary layer properties are presented in § 3.1, and the development of the IBL over the wind farm and the flow-adjustment are discussed in § 3.2.

Boundary layer properties
An overview of the surface forcings and the basic boundary layer properties such as the boundary layer height, friction velocity, jet velocity, and the stability parameter z i /L in the precursor domain, are presented in table 1. We determine the boundary layer height by the method used by Kosović & Curry (2000) and Beare et al. (2006). The boundary layer height z i is defined as the height where the mean stress is 5% of its surface value (z 0.05 ) followed by a linear extrapolation, i.e., z i = z 0.05 /0.95. At higher cooling rates, the friction velocity decreases, which indicates that there is a reduced turbulence in the boundary layer. Furthermore, the boundary layer becomes shallower, i.e. z i reduces.   u 2 + v 2 , where represents planar averaging and the overbar represents time averaging. The strength of the jet, which is defined as the ratio of wind magnitude of the jet to the geostrophic velocity, i.e., u jet /G, increases as the cooling rate increases while the jet height z jet /z h decreases. The jet plays an important role in sustaining continuous turbulence in the boundary layer (Banta 2008;Mahrt 1998). For stronger stratification cases SBL-4 and SBL-5 (see table 1), the height of the jet relative to the turbine hub height is approximately 1.5, which means that the jet height is equal to the height of the top of the turbine blades. Figure 3(b) shows the planar averaged potential temperature profile. The height of the capping inversion z c is defined as the height where the temperature gradient is maximum. For z c /z h 1.8, we observe that the inversion height is approximately equal to the SBL height, such that the direct interaction with the IBL developed by the wind farm is possible. Figure 3(c) presents the wind angle variation α as a function of height for the different cases. For higher cooling rates, a wind veer as strong as 15 • − 20 • is observed.  showed that such a sharp wind veer causes wind turbine wakes to develop a skewed spatial structure.
Based on z i /L, Holtslag & Nieuwstadt (1986) identified three SBLs regimes, namely 1) near-neutral regime (0 < z i /L 1) with weak stability characterized by continuous turbulence, 2) an intermediate regime (1 < z i /L 10) with moderate stability where the boundary layer follows z-less scaling with continuous turbulence, and 3) a highly stable intermittency regime (z i /L > 10) where the turbulence is weak and sporadic and therefore not continuous in time and space. In all the cases considered in the present study, z i /L < 3, indicating weak to moderate stability of the boundary layers. The boundary layer under such conditions remains continuously turbulent, and the similarity theory applies to the surface layer. Furthermore, continuous turbulence is sustained by the high shear of the LLJs.
In addition to z i /L the effect of capping inversion, which takes into account the free atmospheric stratification, can also be characterized by the gradient Richardson number (Ri) calculated by the Brunt-Väisälä frequency N and mechanical shear S: (3.1) where u w = uw + τ xz − u w, and v w = vw + τ yz − v w. The fluxes are normalized with the surface flux of the SBL-1 case to show the reduction in the turbulent momentum flux at higher cooling rates. It is evident from figure 3(d) that the turbulence in the boundary layer reduces when the surface cooling rate is increased. The inset in figure  3(d) shows that the Richardson number Ri increases monotonically with height for all the cases. At the top 10-20% of the boundary layer, the Ri increases above the critical Ri c (based on the hydrodynamic instability theory, Richardson (1920); Taylor (1931); Miles (1986); Galperin et al. (2007)). Zilitinkevich et al. (2008) classify the boundary layer into three regimes: 1) weakly stable regime at Ri < 0.1, 2) a transitional regime at 0.1 Ri 1 with strong turbulence at Ri << 1; and 3) weak turbulence regime at Ri > 1, capable of transporting momentum but not heat. At higher cooling rates (cases SBL-4 and SBL-5), the Ri number increases rapidly with height, limiting the turbulence to very low heights, which affects the IBL dynamics in the presence of a wind farm.
To conclude, the initialization stage yields completely turbulent, quasi-steady boundary layer which serves as an realistic inflow condition for the wind farm.

Flow adjustment in and around the wind farm
Instantaneous flow structures of the horizontal velocity for the SBL-3 case are shown in figure 4. The top panel shows the velocity contours in an x-z plane (through the 3 rd row, note that only the lowest z/D = 5 is shown). We observe significant wake meandering in the vertical direction and that the boundary layer turbulence is restricted to the region below the capping inversion. The meandering of the wakes in the vertical direction can play a significant role in entraining the high-velocity fluid in the jet. The figure shows that the turbine wakes interact directly with the LLJ, inducing mixing, and reducing the strength of the jet. Besides, we note that the LLJ is reduced in intensity after the 5 th column of the wind farm. Due to strong stratification and low lying inversion, the turbulence is more or less limited to the small region above the wind farm. The middle panel of figure 4 shows a horizontal snapshot of the flow at hub height (x-y plane). We notice significant wake meandering in the lateral direction and significant turbulence in the rear of the wind farm. In contrast with a neutral boundary layer with alternating high and low-speed turbulent streaks, the boundary layer, in this case, is only composed of small scale turbulent structures. Finally, the bottom panel shows a y-z plane at a distance (1D) behind the 6 th turbine column. This figure is interesting as it shows a significant spanwise flow of the fluid with the LLJ impinging on the turbine in the 1 st row on the left. This happens due to the wind veer induced by the Coriolis forces and is discussed in § 5. As a result, the turbines in the 1 st row entrain the high-velocity jet, which increases the power production. Another noteworthy point here is that the figure shows the importance of performing non-periodic, fully-finite simulations with a fringe layer in the spanwise direction. If a spanwise fringe layer is not used for the inlet, the The lines indicate the development of the the capping inversion height, and the lines with markers the IBL height as in panel (a). Note that for SBL-4, and SBL-5 the IBL grows above the capping inversion. turbine in the 1 st row would be operating in the wake of the wind farm, consequently under-predicting the turbine power output.
The turbines extract energy from the incoming flow and thereby create a momentum deficit in the wake. The wakes start interacting with the boundary layer both in the lateral and vertical direction via turbulence, and the momentum deficit spreads in the boundary layer, which in turn entrains air towards the turbines. The region of momentum deficit gives rise to the IBL, above which the boundary layer is undisturbed by the dynamics near the surface, whereas inside the IBL, the flow structure changes downwind due to momentum extraction by the turbines. The growth of the IBL shows how the wind farm modifies the flow. Furthermore, the height of the IBL is useful in analytical modeling of the wind farm power production (Meneveau 2012). There is no set rule for calculating the IBL height, for example, Wu & Porté-Agel (2013) define it as the height where the time-averaged wake velocity is 99% of the mean flow velocity at that height, Allaerts & Meyers (2017) define it as the height where the ratio of time-averaged horizontal velocity magnitude and the inflow velocity at the same height, taken in a plane 2 km upwind, reaches a threshold of 97%, Stevens (2016) define it as the height where the vertical energy flux reaches the free stream value. Following Allaerts & Meyers (2017), we define the IBL as the height where the time-averaged horizontal velocity magnitude u mag is 97% of the planar averaged inflow velocity at the same height. Besides, we fix the turbine top (z h +D/2) as the minimum height of the IBL as the IBL grows over the turbine top. Figure  5 (a) shows that the IBL height decreases when the surface cooling rate is increased and increases with the downstream location in the wind farm. The capping inversion height mostly determines the growth of the wind farm IBL. This is analogous to the growth of an IBL over a roughness change due to horizontal advection of air. Here, the wind farm is felt by the upstream flow as a roughness change, and due to the continuity constraint, the flow accelerates over the wind farm.
In a SBL, the capping inversion gets pushed up by the growing IBL (Allaerts & Meyers 2017). In figure 5(b), the base of the capping inversion z c , defined as the height where the temperature gradient is maximum, is plotted along with the IBL for different cases. It is evident from the figure that the capping inversion is pushed up due to the IBL. For the first two cases, the IBL stays below the capping inversion. The displacement of the capping inversion increases with the increased cooling rate, and for SBL-4, and SBL-5, the IBL grows above the capping inversion. The capping inversion is strongly stratified in these cases, and the Ri number of the flow is high at the top of the boundary layer. In these cases, the capping inversion acts as a lid, which limits the growth of the IBL. Due to the continuity constraint, the wind goes around the wind farm. The space between the top of the turbines and the capping inversion height determines how much wind flows around the wind farm. The capping inversion is at the height of z c /z h 2.114 for cases SBL-3, SBL-4, SBL-5, which is approximately 0.5D above the tip of the turbines. In such a scenario, the fluid has very little room to go above the wind farm due to the stabilizing effect of the capping inversion. Consequently, in these cases, we see an appreciable amount of flow going around the wind farm. In essence, the so-called blockage due to the wind farm is the highest for SBL-5 and lowest for SBL-1, which is related to the space between the capping inversion and IBL. Figures 6(a,b) show the time-averaged streamlines at hub height for the cases SBL-1 and SBL-5. Figure 6(a) shows the streamlines for SBL-1; we see that the streamlines are nearly parallel and show marginal divergence. Figure 6(b) shows the streamlines for the SBL-5 case; we observe significant streamline divergence proving that the flow goes around the wind farm. Rominger & Nepf (2011) observe that when a flow encounters the leading edge of a canopy, a part of the flow is diverted, and the remaining part advects through the porous canopy. As the turbines start extracting energy, the shear in the IBL reduces, causing an increase in the Ri number in the IBL. The inset of figure  3(d) shows that the increase in Ri with height is maximum for SBL-5. As the shear in the flow decreases due to the energy extraction by the turbines, the Ri increases. With the increase in local Ri, the flow stability increases, and the fluid finds it difficult to go above the wind farm, and it takes the path of least flow resistance, i.e., around the wind farm. The effect is similar to the flow going around a three-dimensional obstacle like a mountain under highly stratified conditions (Hunt & Snyder 1980;Baines 1979). Figures 6(c,d) show the pressure perturbation normalized by the inlet pressure at the base of the capping inversion z c and at the hub height for the different cases. For SBL-5, the pressure perturbation starts increasing in the entrance of the wind farm when the IBL is at the same height as the capping inversion. As the capping inversion poses resistance to the developing IBL, the flow experiences an adverse pressure gradient. This makes it difficult for the flow to go through or above the wind farm, forcing it to go around.

Energy budget analysis
In the boundary layer, the wind turbines extract energy from the flow and entrain fresh momentum from the upper layers of the atmosphere. An energy budget analysis is a convenient way to understand the diverse phenomena involved in the power production of a wind farm. In § 4.1, a budget analysis of the total energy and its different components is presented, and the turbine power production is discussed in § 4.2.

Entrainment, streamwise flow work
The steady-state, time-averaged, energy equation is obtained by operating the momentum equation with u i (Sagaut 2006 where, the overline represents time averaging, and u i u j = u i u j + τ ij − u i u j represents the momentum flux to which the subgrid scale components have been added. We are interested in the total power production per wind turbine column and energy balance around each turbine. To calculate the total energy, we numerically integrate the terms in (4.1) in a control volume ∀ surrounding each turbine column. Figure 7 schematically represents the dimensions and the extent of the aforementioned control volume. The control volume covers all the turbines in a column and has a streamwise extent of s x D, i.e. 7D, with 3.5D in front and 3.5D behind the turbines, in the streamwise direction. In the vertical direction, the control volume has a dimension of D, and covers the volume 7D 20D D 2.5D 2.5D between z h − D/2 and z h + D/2. In the spanwise direction, the control volume covers the whole column with an additional 2.5D on the sides, essentially 20D. So the total control volume size for each column is, 7D × 20D × D. Integrating (4.1) and rearranging gives P, Turbine power (4.2) In equation 4.2, E k represents the divergence of the kinetic energy flux, which includes both the resolved and the subgrid scale kinetic energy, the turbulent transport term T t , which includes the entrainment of mean momentum due to turbulence and the entrainment of turbulent kinetic energy due to fluctuating velocities (third order terms), T sgs represents the transport of momentum due to SGS fluxes. The flow work F represents the energy transfer due to static pressure drop of the flow across a turbine. The term B represents the turbulence destruction due to buoyancy, G represents the mean geostrophic forcing, and finally P represents the turbine power production. We are interested in the contribution of different budget components to power production. Therefore all the terms are normalized by the magnitude of the power produced by the turbine column. The SGS transport T sgs and the buoyancy fluxes B are small, less than 10% of the first-column power and have been left out of the plots for brevity. The terms in equation (4.2), which include the gradients i.e. E k , T t , T sgs , and F, represent the net flux out of the control volume, for example, E k = E out − E in . Positive values of these terms E in > E out indicate that more energy is added to the control volume than removed. This indicates that in the control volume, energy is extracted from the flow by the turbines or other means. Negative values of these terms indicate E out > E in , which means energy is being added to the flow.
For all the cases, the geostrophic forcing term G remains nearly constant for all the columns of the wind farm, representing a constant driving force. Besides G, there are three primary energy sources, which decide the turbine power production, namely (i) the kinetic energy flux E k , (ii) the work done due to the static pressure drop F, and (iii) the turbulent transport T t , which includes both entrainment of mean momentum into the wind farm by turbulent fluxes (shear production term) and the entrainment due to turbulent fluxes (third-order turbulence terms). Major energy sinks are the power extracted by the turbines P, the dissipation D, and the turbulence destruction due to buoyancy B. Figure 8(a) shows different energy components for the SBL-1 case. The turbines continuously extract energy from the flow, and the kinetic energy flux decreases in the downstream direction. For the last three columns, E k < 0, which means more energy exits the control volume than entering it. This happens because of the entrainment of the kinetic energy T t from fluid above the wind farm. The turbulent transport term T t is composed of fluxes like u·u w and v·v w , which represent the vertical (downward) flux of the mean momentum created by turbulence, i.e., entrainment of mean energy from above towards the turbines. With more and more turbulence created by the turbine wakes, this entrainment flux increases in the downstream direction. In a wind farm operating under neutral stratification and no capping inversion or LLJ, this entrainment flux is of the same order of magnitude as the turbine power production. This flux acts as the major source of power for the downstream wind turbines and reaches a constant value towards the end of the wind farm (Calaf et al. 2010;Cal et al. 2010). A similar variation of energy fluxes has been reported by Allaerts & Meyers (2017) in the simulations of conventionally neutral boundary layers. For SBL-1, the jet height (z jet /z h = 2.670) is well above the wind farm. The IBL grows above the wind farm and facilitates the interaction with the high-velocity jet. Consequently, the entrainment continuously increases downstream and reaches its maximum towards the end of the wind farm. Figure 9 shows that although the jet strength reduces for SBL-1, the jet more or less persists above the entire wind farm. Figure 8(a) shows that the pressure-velocity correlation due to the static pressure drop, also known as the flow work F, is positive and increases along the length of the wind farm. This is an indication that the turbines operate in a favorable pressure gradient in the SBL-1 case. F has a significant contribution towards the power production near the end of the wind farm. The turbine power production, which is the major sink, is maximum at the entrance and reduces downwind due to the effect of the upstream turbine wakes. This variation is typical for a wind farm with an aligned layout and has been observed both in observations and numerical studies (Hansen et al. 2012;Stevens et al. 2014a;Wu & Porté-Agel 2013). The dissipation D acts as an additional energy sink and remains roughly constant as a function of the downstream position in the wind farm. Figure 8(b) presents the energy budget for the SBL-2 case. The figure shows that the entrainment T t increases until the 7 th column when it saturates. A similar trend is observed for the SBL-3 case in figure 8(c), but then the entrainment already saturates after the 5 th column. For SBL-3, the jet height (z jet /z h = 1.836) is slightly above the turbine tip height. The increase and decrease in entrainment correspond to the positions when the wind farm IBL starts interacting with the LLJ. Figure 9 shows that the jet strength for SBL-3 is significantly reduced after the 5 th column. For SBL-5, the jet is utilized by a couple of columns at the entrance, and the remaining columns have little or no jet left to entrain, therefore T t remains nearly constant for this case after the initial increase.
Figure 10(a) shows the variation of D + B for different cases. Both B and D act as energy sinks in the budget, and the buoyancy flux B is small, less than 8% of the firstcolumn power for all the cases, therefore it is combined with the dissipation to represent the net energy sink. B + D is maximum when the turbines interact with the LLJ. This shows that the turbulence production due to mean shear is maximum when the LLJ is at lower heights. In SBL-5, for which the stability is the highest (see figure 8(d)), T t is nearly equal to D, which means there is no effect of entrainment fluxes on the turbine power production and we see a continuous drop in the kinetic energy flux as well as power production. Under stable stratification, increasing the stability damps out the vertical velocity fluctuations, which results in a reduction of in the downward transport of horizontal momentum towards the surface (see figure 3(d)). This results in a reduction of shear production terms u w ∂u/∂z and v w ∂v/∂z in T t , which causes a reduction of the turbulent kinetic energy. As mentioned before, the absolute value of B is not significant. However, the turbulent fluctuations damped out by the stratification, in turn, affect the momentum flux, which causes the weak turbulence in the SBL (Shah & Bou-Zeid 2014). With the jet utilized by the first few turbine columns in SBL-5, the turbines downstream experience a reduction in shear production as well as mean shear. Consequently, we see a continuous decrease in power production for the turbines further downstream. Monin & Yaglom (1971) describe the Obukhov length as the height below which buoyancy or the thermal effects do not play an important role. In a SBL, for z << |L|, the effects of dynamic factors such as shear dominate. For, z > |L| the thermal effects dominate diminishing turbulence. The Obukhov length for cases SBL-5 and SBL-4 are 48.8 m and 66.0 m, respectively, which is less than the turbine hub height. In these cases, the turbines operate mostly in a buoyancy dominated region with high stability. Therefore, we see minimal shear production and turbulent transport T t in these cases. Here, T t is more or less balanced by B+D (figure 8(d)) and the turbine power production depends completely on non-turbulent phenomena such as the divergence of mean kinetic energy flux and the static pressure drop. With the increased shear associated with LLJ, the turbines in a SBL produce more power compared to the turbines operating in the absence of a LLJ. For cases with high stability i.e. z h < |L|, E k , F, and G are the only energy sources available, as T t is balanced by B and D. Therefore, the power production decreases with increasing stratification. However, even in such cases, in the presence of a LLJ, the front turbine columns may still perform well due to the elevated shear in the LLJ. Figure 10(b) presents the variation of the flow work F for different cases. SBL-1, SBL-2, and SBL-3 show that the flow work is always positive, which shows that the turbines operate under a favorable pressure gradient. Since F > 0, it acts as an energy source for the turbine power production for the cases SBL-1, SBL-2, and SBL-3. For the case SBL-5, with the increase in streamwise distance, the resistance to the flow created by the capping inversion increases as the IBL grows. This resistance to the flow reaches a maximum at the third turbine column (approximately x/D ≈ 45) when the IBL height is the same as the height of capping inversion, and we see the minimum of F at this point. Following this critical point, the flow starts going around the wind farm, facilitating the easier flow, and the consequent pressure drop across the wind farm increases.

Turbine power production
Figure 11(a) presents the column-averaged power normalized by the average power production of the turbines in the first column of SBL-1. The figure shows that the power production of the first turbine column increases significantly when the surface cooling is increased. The reason is that the average hub height velocity is higher for the cases with stronger stratification, see figure 3 (a). However, the figure shows that the turbine power production towards the end of the wind farm is lower for cases SBL-4 and SBL-5 than for SBL-3. The reason is that the turbulent energy entrainment further downstream in the wind farm is limited for these cases. To study the effect of wake recovery on the performance of downstream turbine columns for the different cases, figure 11(b) presents the column-averaged power normalized by the power production of the first column. An increase in power production after the second column indicates high turbulence activity and faster relative wake recovery, and a continuous decrease in power indicates slower wake recovery. For SBL-1 the P/P column=1 increases downstream of the first turbine. This increase in the relative power production with the downstream direction indicates that more energy is entrained from the jet which is then extracted by the turbines. For SBL-4 and SBL-5, P/P column=1 decreases asymptotically to a constant value indicating reduced relative wake recovery. Figure 12 (a) and (b) show the turbulent entrainment and wake recovery for different cases. In the region behind the fifth column, the SBL-3 case shows maximum entrainment. For this case, the jet height is z jet /z h = 1.836 and the vertical meandering of the turbine wakes can entrain the high velocity jet, this interaction reaches a maximum around the 3 rd turbine column after which the jet is completely used up and the entrainment continuously decreases. Figure 12(b) shows that SBL-1 has the fastest wake recovery of all the cases. The inlet turbulence intensity at the hub height for this case is the highest with TI z h = 5.82%. Cases SBL-2 and SBL-3 show significant wake recovery towards the end of the wind farm. For these cases the inlet Obukhov length is 189 and 100 meters, respectively, which is greater than the hub height i.e. the turbines are in a regime where there the shear generated turbulence effects dominate, with more and more turbulence generated towards the end of the wind farm, these cases show significant wake recovery. Figure 12(b) shows a significant reduction in the upstream wind velocity in front of the first turbine column, which indicates the effect of the adverse pressure gradient created by the wind farm blockage. This upstream reduction in wind speed increases with stratification and is maximum for SBL-5 for which the adverse pressure gradient caused by the capping inversion is maximum. This flow blockage reduces the inlet wind velocity for the first column of turbines, and the turbines produce lesser power than what they would if they were free-standing. Segalini & Dahlberg (2019), Bleeg et al. (2018), Allaerts & Meyers (2018) and Wu & Porté-Agel (2017) have also reported similar upwind flow reduction due to wind farm blockage.

Effect of wind veer
In the presence of the Coriolis force, the wind follows an Ekman spiral in which the wind velocity vector changes its direction with height. The changes in the wind angle are caused by the imbalance between the pressure gradient and frictional forces. Under stable stratification, the Ekman spiral is very pronounced, and it produces a wind veer facilitated by the significant spanwise flow. In our simulations, we use a PI controller to fix the wind angle at the hub height to zero. This results in a flow that has a positive spanwise velocity below the turbine hub and a negative spanwise velocity above the turbine hub. Figure 13 presents the power map for the SBL-3 case with all the entries normalized by the power produced by the turbines in the 1 st row of their respective columns. It is evident from the figure that the turbines in the 1 st row produce more power compared to the rest of the turbine rows. Furthermore, there is a gradual reduction in power production towards the fourth row. This variation in power is because of the wind veer created by the Coriolis force. We find that this effect is substantial for SBL-3, SBL-4, and SBL-5. The effect is certainly present for SBL-1 and SBL-2 but not significant. Figure 14(a) shows the horizontal velocity magnitude for the SBL-3 case in the y-z plane cut through the middle of the 6 th turbine column. In this case, the jet height is z jet /z h = 1.836, which is slightly above the turbines. We observe that the turbines completely utilize the jet above the wind farm due to entrainment and wake meandering, whereas the jet to the left of the 1 st turbine row provides a continuous supply of fresh momentum due to the spanwise flow which goes to the right. The turbines in the outer rows receive a constant supply of high-speed jet, which is utilized by the turbines, and the remaining fluid goes to the turbines on the right. As the 1 st row has already utilized the jet, the power production of the next row is reduced. Furthermore, the local variation in the wind velocity created by the turbine wakes also causes the wind to deflect clockwise. van der Laan & Sørensen (2017) and Gadde & Stevens (2019) report deflection of the turbine wakes clockwise in the Northern hemisphere due to the imbalance created by the entrainment fluxes induced by the wind farm. We observe a similar clockwise deflection of the turbine wakes as a result of which the turbines in the inner rows operate in the wake of the outer rows. Figure 14(b) shows the streamwise downward energy flux u · u w for the SBL-3 case. The wake structure is skewed due to the lateral shear created by the spanwise flow. The turbine in the 1 st row entrains most energy from the jet, and the subsequent rows entrain less energy from the jet due to the wind turbine wake. This skewed spatial structure of energy entrainment is an additional reason for the observed power variation.

Conclusions
Large-eddy simulations of wind farms in stable boundary layers were carried out in the present study. The objective of the study was two-fold: 1) to study the effect of stable stratification and low-lying capping inversion on the flow development in wind farms, and 2) to study the impact of the LLJs on power production. The study was carried out by systematically varying the cooling rate at the surface, which gives rise to boundary layers with capping inversions and LLJs at various heights.
Decreasing the cooling rate reduces the turbulence near the surface, which in turn increases the geostrophic components in the boundary layer producing LLJs with increased strength. With increasing stratification, the height of the capping inversion reduces, offering increased resistance to the flow. Due to the adverse pressure gradient created by the capping inversion and the blockage created by the presence of a wind farm, the flow goes around the wind farm. When the capping inversion is significantly above the IBL height, the flow accelerates above the wind farm. However, when the IBL height and the capping inversion height are similar, there is significant spanwise flow acceleration around the wind farm. Therefore, performing simulations with periodic boundary conditions in the spanwise direction over-predicts the flow blockage as the flow cannot go around the wind farm. LES of conventionally neutral wind farms by Allaerts & Meyers (2017) and Wu & Porté-Agel (2017) show such flow blockage for 'infinite' wind farms. Furthermore, Allaerts & Meyers (2019) report that the wind farm flow blockage is maximum at an aspect ratio of 3 : 2 (spanwise width-to-streamwise length ratio of the wind farm).
The presence of a LLJ can be beneficial to wind farm power production. When the jet is above the wind farm, it serves as a source of high momentum for the turbines to entrain. We observe elevated entrainment when the jet is above the wind farm. The entrainment is maximum when the turbine wakes can directly interact with the jet by the vertical meandering of the turbine wakes. The entrainment reduces once the turbines entrain all the energy in the jet. If the LLJs are at a height z jet z h + D/2, the outer turbines completely extract the momentum in the LLJ and perform significantly better than the inner turbines. Under similar stability conditions, i.e., stable stratification, a wind farm would perform better if the LLJ is present above the wind farm than when a LLJ is absent. The simulations prove that the outer turbine rows utilize the LLJ and the entrainment shows a downward trend after the jet is completely exhausted. In sites where LLJs are prominent, wind farms with higher aspect ratios (spanwise width-to-streamwise length ratio of the wind farm) are more beneficial than long wind farms with low aspect ratios.
The stable boundary generally has low turbulence intensities, and the surface Obukhov length can serve as an important length scale to predict the impact of the stability. We find that for z h >> |L|, the shear effects dominate, and the entrainment is more than the dissipation and buoyancy destruction. If z h < |L|, the thermal effects dominate, and there is very little entrainment as buoyancy damps out the vertical velocity fluctuations reducing both vertical kinetic energy and downward turbulent fluxes.
In the presence of a LLJ, an appreciable spanwise flow is created by the wind veer. Consequently, the outer turbines produce more power than the inner rows. Inner turbine rows can only interact with the LLJ by the increased turbulent entrainment, while the outer turbine rows can directly extract energy from the LLJ. This effect is prominent only if the jet height z jet ≈ z h + D. Finally, the present study only focuses on the cases where the jet is above the turbine top height i.e. z jet z h + D/2, consequently the turbines only experience positive shear in the LLJ. Further studies are required to analyze the effect of negative shear of the LLJ (when z jet < z h + D/2) on the wind farm power production.

Acknowledgments
We would like to acknowledge Dries Allaerts for the insights into the implementation of the wind angle controller. 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. This work was carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF corporation, the collaborative ICT organization for Dutch education and research. We also acknowledge PRACE for awarding us access to JUWELS based in Germany at Jülich under PRACE project number 2017174146.

Declaration of interests
The authors report no conflict of interest.