Simulation and modeling of wind farms in baroclinic atmospheric boundary layers

As the need for renewable energy sources grows, the wind energy sector resorts to increasingly larger turbines, often clustered in offshore wind farms. These wind farms are regularly built in relatively close proximity to each other, as the available shallow water area suitable for fixed-foundation turbines is limited.^1,2 Consequently, wake losses, both within wind farms^3,4 and between neighboring farms,^1,5,6 are becoming more prominent, severely reducing wind farm power output. It is, therefore, critical to obtain a good understanding of how the interaction between wakes and the atmospheric boundary layer (ABL) governs (large-scale) wake recovery. How this process is affected by various atmospheric phenomena is not yet fully understood and is often unaccounted for in modern wind farm design tools.^7,8

One atmospheric condition whose impact on wind farms remains relatively unexplored is baroclinicity.⁹ In a baroclinic ABL, a horizontal temperature gradient exists throughout the depth of the ABL. As a result, the gradients of pressure and density are not aligned, causing the driving mean pressure gradient (and by extension the geostrophic wind) to vary with height. This contrasts with the barotropic conditions typically assumed when modeling wind farms, in which the geostrophic wind is constant with height. Baroclinic ABLs regularly occur in reality,^10,11 especially at middle latitudes. Typically, baroclinic conditions are found near land-sea transitions,^12,13 water-ice boundaries,¹⁰ around sloping terrain,¹⁴ or near mountain ranges.¹⁵ In fact, even the mean North–South temperature gradient resulting from the uneven heating of the Earth induces a baroclinic effect.¹⁶ 

Field measurements^13,17,18 and numerical simulations^{9,10,16,19–22} indicate that baroclinicity can alter multiple atmospheric properties governing wake–atmosphere interactions. This suggests that baroclinicity can explain various observed discrepancies between measurements and classical engineering wake models. For instance, baroclinicity is known to alter low-level jets,^9,13,15 turbulence levels throughout the ABL,^10,15 wind veer profiles,¹⁰ and boundary layer height,²¹ among other properties.

In previous work,⁹ we employed large eddy simulations (LES) to study how the performance of a wind farm is impacted by a specific type of baroclinicity, called negative geostrophic shear. We showed that geostrophic shear induces momentum fluxes aloft in the ABL and alters the strength and height of low-level jets. This results in reduced downward momentum transport into the wind farm and decreases the relative power output of downstream turbines. However, different forms of baroclinicity yield distinct ABL properties and are expected to affect atmosphere–wind farm interaction differently. Furthermore, the impact of baroclinicity on wind farm wake recovery remains unexplored.

In this work, we build upon the methods and findings presented in Ref. 9 to fill in these knowledge gaps. We evaluate how different types of baroclinicity affect the power production of a wind farm and the wake recovery behind it, as illustrated in Fig. 1. We accomplish this using LES, which offer precise control over atmospheric conditions and enable the isolation of specific effects. LES have been employed extensively in recent years to study atmosphere–wind farm interactions. While initial wind farm LES studies were highly idealized,²³ there has been a shift toward more realistic simulations by incorporation of effects such as thermal stratification,^24–31 Coriolis force,^32–34 dynamic wind direction changes,³⁵ geostrophic shear,⁹ and coupling to mesoscale models.³⁶ 

Furthermore, we introduce a new engineering wake model, developed using presented LES data, designed to predict the large-scale wake recovery behind wind farms. The proposed model provides accurate wind farm wake recovery predictions by incorporating insights into the atmospheric conditions. The model considers wake expansion based on the atmospheric turbulence intensity. Furthermore, it predicts the vertical distribution of the wake deficit by considering the prevailing wind shear, which we identified as crucial for making accurate wind farm wake recovery predictions.

This article is organized as follows. Section II describes the employed LES framework and gives an overview of the considered cases. The boundary layer characteristics are discussed in Sec. III. Section IV discusses the flow in and around the wind farm. In Sec. V, the wind farm performance is discussed, and in Sec. VI, we evaluate the wind farm wake recovery. The newly developed engineering wind farm wake model is introduced in Sec. VII, and comparisons are made with the LES results. Finally, the conclusions are presented in Sec. VIII.

Note that Eq. (3) is solved for a modified potential temperature $θ̃∗$ and contains an additional advection term $−ũΘ|x−ṽΘ|y$ compared to traditional LES, which are unique to baroclinic simulations. Additionally, evaluation of the temperature perturbation $θ̃′$ in Eq. (2) requires special attention in baroclinic conditions. These terms are separately discussed in Sec. II B.

The total subgrid-scale (SGS) stress tensor is decomposed as $τijt=τij+(τkkt/3)δij$⁠, where $τij=uiuj̃−ũiũj$ denotes its deviatoric part, and the trace $τkkt/3$ is absorbed into the modified pressure, defined as $p̃∗=(p̃−p∞)/ρ0+τkkt/3$⁠, with mean pressure $p∞$ and hydrostatic density $ρ0$⁠. $δij$ is the Kronecker delta. Stresses due to molecular viscosity are neglected as is customary for high Reynolds number LES, and $qi=uiθ̃−ũiθ̃$ is the SGS heat flux vector. $τij$ and $qi$ are modeled using the anisotropic minimum dissipation (AMD) model,^39,40 which was found suitable for modeling heterogeneous turbulent flows with relatively limited computational overhead.³⁸ 

The geostrophic wind velocity is defined as $Ui=−ϵij3∂jp∞/(ρ0fc)$⁠, with alternating unit tensor $ϵij3$ and Coriolis frequency $fc$⁠. Finally, the turbine force $fi$ is modeled using the actuator disk model,^23,41,42 as described in  Appendix A, which is well validated and allows large-scale interactions between wind farm and ABL to be captured accurately.^43,44

The LES code is pseudo-spectral with periodic horizontal boundary conditions. Consequently, we cannot directly include a global horizontal temperature gradient, as that would violate the horizontal periodicity of the resolved potential temperature field $θ̃$⁠. Therefore, $θ̃$ is decomposed into two components as $θ̃(x,y,z,t)=θ̃∗(x,y,z,t)+Θ(x,y)$⁠. Here, $Θ(x,y)=Θ|xx+Θ|yy$ denotes the baroclinic temperature contribution, assuming that $Θ|x$ and $Θ|y$ are constant in space and time (and by extent $U|z$ and $V|z$⁠). We solve for the modified potential temperature $θ̃∗$ instead of $θ̃$ in Eq. (3), which does satisfy the periodic boundary conditions.

The temperature perturbation $θ̃′$ in Eq. (2) cannot easily be evaluated directly in baroclinic conditions, as $θ̃$ can vary globally in both space and time. However, per our problem definition, global planar heterogeneity in $θ̃$ arises only due to the baroclinic temperature contribution $Θ$⁠. Thus, we may evaluate the temperature perturbation as $θ̃′=θ̃*−⟨θ̃*⟩$⁠, where $⟨.⟩$ denotes a planar average.

We employ a free-slip boundary with zero vertical velocity at the top of the domain. At the wall, we use the Monin–Obukhov similarity theory⁴⁵ to model the stresses, see  Appendix A. The heat flux at the bottom of the domain is set to zero, corresponding to neutrally stratified conditions.

The numerical domain is discretized using a rectilinear grid, which is uniform in streamwise and spanwise directions and comprises $nx=Lx/Δx$ and $ny=Ly/Δy$ grid points, respectively. $Lx$ and $Ly$ denote the domain length and width and $Δx$ and $Δy$ are the grid resolutions. In the vertical direction, the grid is uniform up to a height $zu$ with a grid spacing $Δz$⁠, and above it is stretched using a hyperbolic tangent stretching function.

Derivatives in the vertical direction are calculated using a second-order central finite difference scheme, while pseudo-spectral differentiation is used to compute derivatives in the horizontal directions. We employ the 3/2 anti-aliasing method to prevent aliasing errors in the non-linear terms. A third-order accurate Adams–Bashforth scheme is used for time integration.

A concurrent precursor method⁴⁶ is employed to ensure the inflow region is not contaminated by remaining wind farm flow structures, due to the horizontal periodic boundary conditions. Furthermore, to reduce the reflection and the fringe-induced excitation of gravity waves triggered by the wind farm, a Rayleigh damping layer⁴⁷ is applied in the upper part of the domain. A schematic illustration of the computational domain is shown in Fig. 2.

Note that tildes denoting LES filtering and the asterisk denoting modification of variables are omitted in the remainder of this article for convenience.

While S and $β$ span a full 2D parameter space, we explore the effects of baroclinicity considering four baroclinic cases, based on the orientation of the temperature gradient $∇Θ$ (and by extent the geostrophic shear) with respect to the surface geostrophic wind.^16,20–22 In addition, we include a barotropic reference case. Choosing $Ui(z=0)=(U0,0,0)$⁠, the cases are

• Barotropic (⁠ $S=0, β$⁠),

• Positive geostrophic shear (⁠ $S>0, β=0°$⁠),

• Cold advection (⁠ $S>0, β=90°$⁠),

• Negative geostrophic shear (⁠ $S>0, β=180°$⁠),

• Warm advection (⁠ $S>0, β=270°$⁠).

In the barotropic case, there is no horizontal temperature gradient and the geostrophic wind is constant with height. For positive and negative geostrophic shear, the magnitude of $Ui$ varies with height while its direction remains unchanged. In these cases, $∇Θ$ is perpendicular to $Ui$⁠, such that there is no geostrophic heat advection, see Eq. (5). For warm and cold advection, the geostrophic shear is perpendicular to the surface geostrophic wind $Ui(z=0)$⁠; thus, $Ui$ changes both in magnitude and direction with height. Over time, the simulation domain experiences heating or cooling due to geostrophic heat advection. This ensues from a non-zero geostrophic wind component in the direction of the temperature gradient $∇Θ$⁠.

For all baroclinic cases, we set the geostrophic shear magnitude to $6 ms−1km−1$⁠. Similar values have been observed in field measurements.^13,49 For example, Floors et al.¹³ found that notable thermal winds were present in the data from a measurement campaign in Høvsøre, at the Danish coast. The geostrophic shear component related to the land-sea transition showed values of up to $10 ms−1km−1$ in extreme cases, while the mean value and standard deviation appeared to be roughly 2 and $3 ms−1km−1$⁠, respectively.

An overview of the considered cases is given in Table I. We use a reference temperature of $θr=286$ K and a Coriolis frequency of $fc=1.159×10−4 s−1$⁠, corresponding to North Sea conditions, meaning the considered shear corresponds to a temperature gradient of $|∇Θ|≈2 K$ per $100 km$⁠. Figure 3 presents the cases in the baroclinic parameter space.

The surface geostrophic wind $Ui(z=0)$ varies per simulation (see Table II) and is set to ensure that the hub height velocity is $uh≈9 m/s$ for all cases. Due to the Ekman spiral, the wind direction is typically not directed in the x direction. Therefore, the computational domain is defined with respect to the geostrophic wind vector $Ui$ such that the velocity at hub height is approximately aligned with the x direction.

We consider a numerical domain of streamwise length $Lx=102.4 km$⁠, spanwise width $Ly=10.24 km$ and height $Lz=10 km$⁠, which is discretized using $2048×512×384$ grid points. The resulting streamwise and spanwise grid resolutions are $Δx=50 m$ and $Δy=20 m$⁠, respectively. Up to $zu=1.5 km$⁠, a resolution of $Δz=10 m$ is employed, while above the grid is stretched up to a maximum spacing of $Δz=62 m$⁠. The fringe region for the concurrent precursor method spans the latter $3 km$ of the domain in streamwise direction and $1 km$ in spanwise direction. The precursor simulation and successor simulation are run on the same numerical grid. Both streamwise and spanwise fringes are used, since the Ekman spiral can cause large wakes to significantly deflect in the lateral direction.

The surface roughness height is set to $z0=0.002 m$⁠, typical for North Sea conditions.²⁸ The initial temperature profile is a neutral layer with $θ0=286 K$ up to $z=1 km$⁠, capped by an inversion of $3 K$ over $200 m$⁠, and a free atmosphere stratification of $5 K/km$ (see Fig. 4). Random perturbations are added in the lower portion of the ABL to spin up turbulence.

In all simulations, we consider a wind farm with ten rows and six columns of turbines with an aligned layout, see Fig. 2. The streamwise and spanwise turbine spacing are set to, respectively, $sx=7D$ and $sy=5D$⁠. The turbine diameter and hub height are set to $D=178 m$ and $zh=119 m$⁠, modeled after the DTU 10 MW reference turbine.⁵⁰ The thrust coefficient is set to $CT=0.75$⁠, and the induction factor to $a=0.25$⁠. In the simulations, the actuator disks are automatically yawed to be perpendicular to the local incoming wind.

All simulations are performed in two stages. First, a spin-up simulation is performed on a coarser domain discretized by $nx/2×ny/2×nz$ grid points. After $7 h$⁠, once the ABL reaches a quasi-steady state, the flow field is interpolated to the full resolution of $nx×ny×nz$⁠, and the turbines are inserted. Transient effects from the initial condition were found to be diminished at this point. After an $1 h$ adjustment period, statistics are collected from the 8th to the 11th hour. More details on the time-averaging window are given in  Appendix B.

Figure 4 displays the ABL properties as a function of height. The mean horizontal wind magnitude $uABL=⟨u¯2+v¯2⟩$ is given in Fig. 4(a). Here $.¯$ indicates a temporal average. In the surface layer, the wind velocity is primarily governed by friction with Earth's surface. Consequently, despite the strong variation in geostrophic forcing, the velocity profiles in the surface layer overlap for all cases. However, at higher elevations, starting from roughly a third of the ABL height, baroclinicity significantly alters the wind velocity. Whether it increases or decreases depends on the type of baroclinicity. Positive Shear and Warm Advection increase the wind speed in the upper ABL, while Negative Shear decreases it. For Cold Advection, we find that the wind speed profile is not significantly affected. Above the capping inversion, the wind speed approaches the enforced geostrophic wind speed as expected.

Figure 4(b) shows the mean wind angle $α=⟨ tan−1(v¯/u¯)⟩$ with height. Note that all simulations have $α(zh)=0°$ by design. Wind veer is governed by Earth's rotation and typically negative in the Northern Hemisphere, as observed for the Barotropic case. Interestingly, baroclinicity inverts the direction of the wind veer in the Ekman spiral for Positive Shear and Cold Advection. For Negative Shear, the magnitude of the wind veer is increased relative to Barotropic.

Figure 4(c) displays the mean turbulent stresses $τ=⟨(u′¯w′¯)2+(v′¯w′¯)2⟩$⁠, with $u′¯w′¯=u¯w¯+τ¯xz−u¯w¯$ and $v′¯w′¯=v¯w¯+τ¯yz−v¯w¯$⁠. In contrast to the wind speed and wind veer, the turbulent stresses are already impacted significantly by baroclinicity at hub height. The turbulent stress decreases monotonically in the Barotropic case. The baroclinic simulations, except Positive Shear, exhibit a local maximum within the ABL. The stresses decrease to zero at the capping inversion height, above which the flow is mostly laminar.

The profiles of the mean potential temperature $⟨θ¯⟩$⁠, relative to $θ0$⁠, are shown in Fig. 4(d). Since we consider conventionally neutral ABLs, $⟨θ¯⟩$ remains constant up to the capping inversion height. The Cold and Warm Advection cases differ from the Barotropic case as the alignment between temperature gradient and (geostrophic) wind direction effectively cools or heats the boundary layer through thermal advection. This observation aligns with findings from previous studies on the effects of baroclinicity in ABLs.^10,16,21 Note that the temperature variation in time induced by the geostrophic heat advection depends on the averaging window chosen.

Table II shows the geostrophic wind at the surface $|Ui|0$⁠, friction velocity $u∗$⁠, and horizontal turbulence intensity $IABL=⟨u′2¯+v′2¯/u¯2+v¯2⟩$ at hub height. $|Ui|0$ differs per case to ensure that the hub height velocity is equal for all cases. The friction velocity $u∗≈0.33 m/s$ is approximately constant, since the wind speed profile in the surface layer is mostly unaffected by baroclinicity. The hub height turbulence intensity varies similarly to the shear stresses, increasing in Positive Shear and Warm Advection conditions and decreasing in Negative Shear and Cold Advection conditions, compared to the Barotropic case.

Baroclinicity affects the flow in and around a wind farm in multiple ways. Figure 5 shows the instantaneous velocity magnitude at hub height for the Barotropic, Positive Shear, and Negative Shear cases. Differences in incoming turbulent structures can be observed for the three simulations. For Positive Shear, larger turbulent structures are visible, while turbulent structures are smaller under Negative Shear. The formation of turbine wakes within the wind farm is also clearly observable. The wake behind the wind farm is affected significantly by baroclinicity. It recovers quicker for Positive Shear, relative to the Barotropic case. Under Negative Shear, the wake recovers slower and exhibits a much stronger wake deflection. This is expected, as wind farm wakes recover primarily due to downward entrainment of kinetic energy from layers higher up in the atmosphere.⁴ For Negative Shear, the undisturbed ABL exhibits increased wind veer with height, see Fig. 4(b), and thus, wake deflection is enhanced.

The variations in the wake recovery rate can be explained by considering the vertical entrainment as well. The entrainment of kinetic energy is significantly affected by the atmospheric turbulence intensity⁴ and velocity resources aloft.^9,51 Figure 6 displays the instantaneous velocity magnitude at the first turbine column at $y=2.3 km$⁠. The contour plots visualize the vertical structure of the atmosphere above the wind farm. Under Positive Shear, the velocity resources and turbulence above the wind farm are higher than in the Barotropic case, see also Fig. 4(a). Consequently, the interaction between the wakes in the wind farm and the atmosphere above is enhanced. The opposite effect is observed for Negative Shear, as the wakes are then mostly confined to the surface layer. Thus, based on atmospheric conditions depicted in Fig. 4, at hub height wakes are anticipated to recover slower under Negative Shear and faster under Positive Shear conditions. How the recovery of the wind farm wake is quantitatively governed by the ABL properties is further examined in Sec. VI.

Figure 7 shows the spanwise-averaged velocity magnitude at hub height $us=⟨u¯2+v¯2⟩$ inside the wind farm. Here, the average is taken over the full spanwise extent of the wind farm. We observe the expected sudden reductions in flow velocity at the turbine locations, followed by gradual recovery behind the turbines. The inflow velocity profile in the surface layer is consistent across all simulations. Consequently, variations in the wind farm entrance region are minimal. However, further downstream the flow through the wind farm is altered in baroclinic conditions. In particular, the wake recovery rate is decreased in negative shear conditions.

Figure 8 shows the streamwise development of the spanwise-averaged horizontal turbulence intensity $⟨I⟩$⁠. In front of the wind farm, $⟨I⟩$ is equal to the value for undisturbed inflow documented in Table II. The first turbine row adds around 4% additional turbulence intensity. Further downstream $⟨I⟩$ approaches an equilibrium state, fluctuating around roughly 20%. Note that the reported turbulence intensity is averaged over the spanwise extent of the wind farm; significantly higher turbulence intensities are found locally. Near the end of the wind farm, the atmospheric turbulence constitutes less than half of the total $⟨I⟩$⁠, and the wake-generated turbulence has partly compensated for the baroclinicity-induced differences in turbulence intensity. Consequently, differences in atmospheric turbulence are expected to be mostly of importance in the wind farm entrance region.

Figure 9(a) shows the column-averaged aerodynamic power $P$⁠, normalized by the first row power $P1$⁠, as a function of turbine row. We find that the overall trend in power production is similar for all simulations. The first row produces the highest power output, benefiting from an undisturbed inflow. The second row sees a dramatic decrease in power production due to wake losses. Subsequently, turbulent transport kicks in and starts replenishing energy taken out by the turbines, via energy entrainment from the surrounding high-velocity wind. Further downstream in the wind farm these processes are roughly in equilibrium with each other and the power approaches a constant value.

To focus on how power production is affected by baroclinicity, we normalize the results with the Barotropic case in Fig. 9(b). The first-row power is roughly the same in all simulations. This is expected, as by design the undisturbed inflow always has a hub height velocity of $uh≈9 m/s$⁠. The production of the second row is strongly impacted by atmospheric baroclinicity. For Warm Advection, this row produces 20% more power, while for Cold Advection there is a 13% power reduction, compared to barotropic conditions. For Positive and Negative Shear, a minor power increase is observed. In addition to this entrance effect [in region I in Fig. 9(b)], a secondary effect exists further downstream (i.e., in region II), changing the trend. For Cold Advection, power production nearly recovers to the Barotropic reference, while for Negative Shear, the power production reduces downstream. For Positive Shear and Warm Advection, notable power increases are observed, of circa 5% and 10%, respectively.

We focus only on the mechanical energy flux $Em=Ek+F$ (we note that the turbulent component to the flow work is negligibly small) and turbulent transport $Tt$⁠, shown in Fig. 10, which are the dominant energy sources to the wind farm. The mechanical energy flux is largest for the first turbine row, since it experiences undisturbed inflow. From the second row onward, it drops significantly as wake losses become more prominent. The development of $Em$ shows little dependence on baroclinicity. This is expected, as the mean inflow profile and flow direction around hub height are highly similar in all cases, as depicted in Fig. 4.

The turbulent transport $Tt$ shows clear variations with baroclinicity. At the wind farm entrance, region I in Fig. 10(b), $Tt$ starts off low, as there is little wake-generated turbulence. As the velocity resources in the surface layer are similar for all cases, the initial differences in $Tt$ result from differences in atmospheric turbulence; see Fig. 4(c). Further downstream, in region II, wake generated turbulence dominates, so differences in atmospheric turbulence are less important. Instead, the availability of velocity resources aloft in the ABL, see Fig. 4(a), governs the turbulent transport $Tt$⁠. The Cold Advection and Barotropic cases exhibit similar velocity resources aloft, and therefore similar $Tt$ downstream. For Warm Advection and Positive Shear, $Tt$ is increased relative to the Barotropic case. For Negative Shear, it is decreased.

In summary, we can explain the effect of baroclinicity on the power output in Fig. 9(b) considering two effects. First, there is an entrance effect (in region I), which results from a combination of differences in atmospheric turbulence, wind veer, and pressure gradients. Second, there is a downstream effect (in region II) governed by the availability of velocity resources aloft in the ABL governing the downward turbulent entrainment of energy.

Baroclinicity significantly affects the recovery of wind farm wakes. Figure 11(a) shows the spanwise-averaged velocity $us$ as a function of distance behind the wind farm $xe=x−xend$⁠. Here, $xend$ is the location “immediately behind” the wind farm, i.e., at a distance of 2D behind the rearmost turbines.⁵⁴ As the wake recovery of large wind farms is mainly governed by entrainment of energy from above, the variations in turbulent entrainment in region II of Fig. 10(b) are also reflected in $us(zh)$⁠. For Warm Advection and Positive Shear, the wake recovery is faster than in Barotropic conditions, while under Negative Shear conditions, it is significantly slower. This is further visualized in Fig. 11(b), which shows hub height contours corresponding to velocity recovery to 95% of the undisturbed flow. For Negative Shear conditions, the wind farm wake length is almost twice as long as in the Barotropic case. In all simulations, the wind farm wake shape in the spanwise direction falls between a Gaussian and a top-hat shape.

Figure 12 displays vertical profiles of the velocity deficit $δ=uABL−us$⁠, for different distances $xe$ behind the wind farm. These vertical deficit profiles exhibit near-Gaussian shapes. The development of the deficit profiles with increasing $xe$ reveals that hub height wake recovery is governed by two mechanisms. First, there is wake expansion. As $xe$ increases, the wake increasingly mixes with the surrounding high-velocity air due to turbulence. Consequently, the initial velocity deficit at $xe=0$ spreads over a larger vertical extent such that the velocity at hub height recovers. Consistent with existing literature,^4,7 simulations with higher atmospheric turbulence intensity (refer to Table II) show faster wake expansion, and the opposite is true for simulations with lower turbulence intensity. Second, we observe an upward displacement of the wake center deficit $zc$⁠; the “bulk” of the wake develops to be situated above hub height. Note that this is not a result of advection of the wind farm wake. Rather, the vertical entrainment (and consequently the wake recovery rate) inside the wake depends on the elevation. Air layers in the upper part of the wake lose energy to lower situated layers at a higher rate than they gain energy from higher layers. As a result, the height $zc$ at which the deficit $δ$ is largest, increases with $xe$⁠. Since the ABL exhibits significant wind shear, the absolute velocity $us$ remains larger at the wake center than below. Consequently, energy remains being entrained downward and away from the wake center to the lower part of the wake. The wind shear in the undisturbed ABL, thus, enables the mechanism driving the wake center upwards. We observe that the rate at which the wake center displaces upwards (Fig. 12) scales with the velocity shear in the ABL, see Fig. 4(a). The effect is significantly less pronounced under Negative Shear conditions, for which velocity resources aloft are lower. Consequently, the center of the wake remains very close to hub height, resulting in a significantly slower hub height wake recovery than for the other cases.

From Fig. 12, it has become apparent that predicting large-scale wind farm wake recovery with a simplified model requires both wake expansion and vertical wake center displacement to be considered. Existing wake models^55–58 are typically designed to predict the velocity deficit inside wind farms. Several of these models consider the effects of atmospheric turbulence on wake expansion.^57–60 However, they do not account for larger-scale interactions with the atmosphere, such as the wake center displacement effect, which makes them less suitable for predicting wind farm wake recovery.⁸ 

Therefore, we propose a simple new model for wind farm wake recovery, in which the wake is considered as one single structure. Based on observations from our LES, we model the velocity deficit in the wind farm wake to follow a Gaussian distribution vertically and to be uniform in the spanwise direction, approximating an extended wind farm, as illustrated in Fig. 13. Since vertical entrainment is the main factor in wind farm wake recovery, spanwise variations within the wake are not considered. Rather our model prioritizes accurate representation of the vertical wake structure. To this end, the velocity and turbulence profiles within the ABL are exploited.

The factor $g(xe)$ accounts for the ground effect by conserving the integrated wake deficit above ground (i.e., for $z≥0$⁠), to ensure conservation of mass. Note that the part of the Gaussian that remains above the ground varies with $xe$ both due to the wake expansion and the displacement.  Appendix C provides more details on the computation of $g(xe)$⁠.

In summary, the proposed model aims to accurately predict the large-scale development of the wind farm wake deficit $δ(xe,z)$ based on the prevailing atmospheric conditions. It exploits insights into the velocity [ $uABL(z)$] and turbulence intensity [ $IABL(z)$] profiles. The model requires the initial velocity deficit of the wake $δ0$ as input, which for the considered simulations is in the range of $0.2−0.25uh$⁠. When not available, $δ0$ can be estimated using conventional models suitable for predicting the flow inside a wind farm. Finally, we note that the employed empirical formulations could be further generalized, which remains a topic for future research.

Figure 11(a) shows a comparison of the recovery of the spanwise-averaged velocity at hub height $us(zh)$ with distance behind the wind farm $xe$⁠, between LES and the proposed model. The model agrees well with the LES and accurately reflects differences in wake recovery rate under different baroclinic conditions, due to consideration of the wind shear and atmospheric turbulence. However, far down stream, we find the model starts underpredicting the wake recovery rate, most clearly observable for the Positive Shear and Warm Advection cases. This discrepancy possibly results from phenomena not accounted for in the model, such as spanwise wake deflection and geostrophic forcing.

As discussed in Sec. VII B, the proposed model aims to accurately predict the vertical structure of the wind farm wake. Figure 12 shows the model-predicted vertical deficit profiles for different distances behind the wind farm. The effects of baroclinicity on both the wake expansion and the wake displacement are well captured by the model. For Positive Shear and Warm Advection, the additional velocity shear and turbulence intensity in the ABL enhance wake recovery relative to Barotropic conditions. For Negative Shear, the velocity shear in the ABL is significantly lower. Consequently, the wake center remains closer to hub height and wake recovery is reduced. For Cold Advection, the wake expansion rate is slower than in Barotropic conditions, as atmospheric turbulence levels are relatively low. However, the bulk of the wake deficit mainly exists aloft in the ABL, such that the difference in wake recovery at hub height is negligible. This highlights the necessity of considering both velocity shear and turbulence intensity in predicting wind farm wake recovery.

For very large distances behind the wind farm, i.e., $xe≥16km$⁠, we find that the model deficit starts to show decreased recovery with respect to the LES, despite not showing a lower wake expansion. This further strengthens the hypothesis that phenomena such as spanwise wake deflection, and additional wake recovery mechanisms, are a source of discrepancies between LES and model.

By employing large eddy simulation (LES), we evaluated the effects of baroclinicity on the performance of a large wind farm and the large-scale wake recovery behind it. Our findings indicate that baroclinicity affects both the velocity shear and atmospheric turbulence in the atmospheric boundary layer (ABL), which in turn impact the wind farm. We considered four baroclinic ABL scenarios: Positive Shear, Negative Shear, Cold Advection, and Warm Advection, and compared them to a Barotropic reference case.

The impact on wind farm performance can be separated into an entrance effect and a downstream effect. Baroclinicity significantly influences the initial wake recovery after the first turbine row. For Warm Advection, the second row experiences a 20% power increase relative to the Barotropic scenario, while for Cold Advection the power production is reduced by 13%. Further downstream, the power production is predominantly governed by turbulent entrainment.

We found that the large-scale wake recovery at hub height behind the wind farm can primarily be attributed to two key mechanisms related to turbulent entrainment. First, the expansion of the wake, governed by the prevailing turbulence intensity levels. Second, an upward displacement of the center of the wake deficit, resulting from height-dependent entrainment of kinetic energy inside the wake, which is enabled by the velocity shear in the ABL. We propose a simple engineering model to predict wind farm wake recovery based on these two mechanisms. The model accurately reflects the LES results for both barotropic and baroclinic scenarios, as it accounts for the relevant atmospheric conditions.

Future research should investigate the representativeness of the model's empirically determined constants under varying atmospheric conditions, such as different thermal stratifications and wind farm configurations, and explore to what extent these constants can be predicted a priori. Additionally, more research is recommended to explore how baroclinicity impacts blockage and wind farm-generated gravity waves,^29,52 which has not been studied in this work.

This project has received funding from the European Research Council under the Horizon Europe program (Grant No. 101124815). We acknowledge the Dutch national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research. We also acknowledge EuroHPC JU for awarding this project access to Discoverer, hosted by Sofia Tech Park. Furthermore, we acknowledge PRACE for awarding us access to IRENE at Très Grand Centre de Calcul (TGCC) du CEA, France (Project 2021250115).

The authors have no conflicts to disclose.

J. H. Kasper: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). A. Stieren: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Writing – review & editing (supporting). R. J. A. M. Stevens: Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Resources (lead); Software (equal); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Due to the inclusion of the Coriolis force in our simulations, the ABL is subject to inertial oscillations, with a period of $2π/fc≈15 h$⁠. Ideally, the time-averaging window would span a full inertial period. However, this was found unfeasible, given the computational resources required. Instead, we have chosen a smaller window that starts well after the initial transient effects have diminished. During this window, the friction velocity remained relatively constant, as shown in Fig. 14, indicating suitability for collecting the statistics.

The empirical formula presented in Eq. (13) was found to work sufficiently well for the cases presented in this work. However, it is not necessarily universal and may need special attention in more extreme conditions. For example, for strong low level jets, which may be found in stable ABLs, the growth of $zc$ might need to be truncated.