Investigating the physical mechanisms that modify wind plant blockage in stable boundary layers

. Wind plants slow down the approaching wind, a phenomenon known as blockage. Wind plant blockage undermines turbine performance for front-row turbines and potentially for turbines deeper into the array. We use large-eddy simulations to characterize blockage upstream of a finite-size wind plant in flat terrain for different atmospheric stability conditions, and investigate the physical mechanisms modifying the flow upstream of the turbines. To examine the influence of atmospheric stability, we compare simulations of two stably stratified boundary layers using the Weather Research and Forecasting model 5 in large-eddy simulation mode, representing wind turbines using the generalized actuator disk approach. For a wind plant, a faster cooling rate at the surface, which produces stronger stably stratified flow in the boundary layer, amplifies blockage. As a novelty, we investigate the physical mechanisms amplifying blockage by evaluating the different terms in the momentum conservation equation within the turbine rotor layer. The velocity deceleration upstream of a wind plant is caused by an adverse pressure gradient and momentum advection out of the turbine rotor layer. The cumulative deceleration of the flow upstream of 10 the front-row turbines sets in motion a secondary circulation. The horizontal flow is diverted vertically, reducing momentum availability in the turbine rotor layer. Although the adverse pressure gradient upstream of the wind plant remains unchanged with atmospheric stability, vertical momentum advection is amplified in the more strongly stable boundary layer, mainly by larger shear of the horizontal velocity, thus increasing the blockage effect.


Introduction
Wind turbines and wind plants each modify the flow within their vicinity.The flow downstream (i.e., the wake) is slower and more turbulent as turbines extract kinetic energy from the wind.The wind also decelerates upstream of a wind plant, an effect known as blockage (Wu and Porté-Agel, 2017;Bleeg et al., 2018;Segalini and Dahlberg, 2020;Sebastiani et al., 2021;Schneemann et al., 2021;Bleeg and Montavon, 2022;Strickland and Stevens, 2022;Jacquet et al., 2022;Strickland et al., 2022).This region of slower wind speed is called the induction region of the wind plant.Wind turbine and wind plant wakes can also affect power production of downstream turbines and plants, an effect know as wake loss (e.g., El-Asha et al., 2017).Wake losses are accounted for in energy-production loss estimates (Filippelli et al., 2018)  annual energy production (AEP) in wind plants (Lee and Fields, 2021).Blockage also reduces power production in wind plants (Bleeg et al., 2018;Sebastiani et al., 2021).However, blockage is usually neglected in loss estimates due to uncertainty in the magnitude of this effect, possibly resulting in lower-than-forecasted energy production and financial losses for wind plant operators (Brower, 2012;Bleeg et al., 2018;Ørsted, 2019;Lee and Fields, 2021).
There is disagreement on the magnitude of the velocity deceleration within the induction region of wind plants.The velocity deceleration within the induction region can vary substantially depending on the size and layout of the wind plant, atmospheric conditions, wind turbine characteristics, and wind speed (Allaerts andMeyers, 2017, 2018;Bleeg et al., 2018;Bleeg and Montavon, 2022;Schneemann et al., 2021;Centurelli et al., 2021;Segalini and Dahlberg, 2020;Strickland et al., 2022).A limited set of simulations report changes in hub-height wind speed larger than 10% at a distance of 2 rotor diameters (2D) upstream of the first row of turbines (Allaerts and Meyers, 2017;Wu and Porté-Agel, 2017;Allaerts and Meyers, 2018).
Conversely, experimental studies and numerical simulations report between 1% and 5% changes in hub-height wind speed at a distance of 1.5D to 3D upstream of the first row of turbines (Bleeg et al., 2018;Bleeg and Montavon, 2022;Schneemann et al., 2021;Centurelli et al., 2021;Segalini and Dahlberg, 2020;Jacquet et al., 2022;Strickland et al., 2022).The large spread in the magnitude of the blockage effect suggests there may be multiple physical mechanisms modifying the velocity deceleration in the induction region.
Idealized simulations of large wind plants suggest gravity waves could potentially amplify blockage (Wu and Porté-Agel, 2017;Allaerts andMeyers, 2017, 2018;Maas, 2022).Large-eddy simulations (LES) by Wu and Porté-Agel (2017) show increasing stratification in the troposphere can result in upstream propagating gravity waves.However, gravity waves do not form in their weak free-atmosphere stratification simulation (Wu and Porté-Agel, 2017).Stronger surface layer stability, longer wind plants, and shallower boundary layers may also increase gravity wave-amplified blockage (Allaerts andMeyers, 2018, 2017;Maas, 2022).Note that Allaerts andMeyers (2017, 2018); Maas (2022) simulate the flow around an infinitely wide wind plant; therefore, the large vertical boundary layer displacement that excites gravity waves (and thus the velocity deceleration in the induction region) is likely overestimated compared to operational wind plants.Wind plants approach the infinitely-wide regime when they are two orders of magnitude wider than the boundary layer height (Allaerts and Meyers, 2019); hence, a 100 km wide wind plant can be considered infinitely wide for a 1 km deep boundary layer.Even though gravity waves in idealized simulations can potentially produce O(10%) velocity reductions in the induction region (x < −2D), these large decelerations have not yet been observed in operational wind plants.
The nature of the physical mechanism modifying blockage in the absence of gravity waves, where the velocity slowdown is on the order of 1-5% of freestream (Bleeg et al., 2018;Segalini and Dahlberg, 2020;Centurelli et al., 2021;Schneemann et al., 2021;Bleeg and Montavon, 2022;Jacquet et al., 2022), has not been studied in depth.Atmospheric conditions have been shown to modify blockage, though the mechanisms through which this happens are unclear.Using a scanning lidar, Schneemann et al. (2021) measured 4% ± 2% velocity reduction between 30D and 5D upstream of the first turbine row during stable surface conditions.They did not observe blockage during unstable surface conditions (Schneemann et al., 2021).
Similarly, Bleeg and Montavon (2022) show larger power reductions for a single turbine row with increasing surface-layer and free-atmosphere stratification.Simulating an infinitely wide wind plant, Strickland et al. (2022)  blockage with stronger surface layer stability.They propose blockage in stable conditions is amplified by a cold air anomaly that produces a high pressure region upstream of the wind plant (Strickland et al., 2022).However, they do not quantify the increase in horizontal pressure gradient caused by the cold air anomaly nor its relative importance compared with other forcing mechanisms (Strickland et al., 2022).Wind plant layout also influences the velocity deceleration in the induction region.
Using neutral LES, Strickland and Stevens (2022) show an increase in the adverse pressure gradient upstream of wind plants with closely-spaced turbines in the cross-stream direction.Bleeg and Montavon (2022) also show blockage varies between an infinitely-wide and finite-sized wind plant.The majority of studies on blockage focus on the parameters (e.g., turbine layout, atmospheric conditions) that modify blockage but not the physical mechanisms that influence the flow.
Here, we investigate how atmospheric stability modifies upstream blockage in the absence of gravity waves.Specifically, we investigate: (1) if the wind speed deceleration upstream of a finite-sized wind plant is modified with atmospheric stability, and (2) the physical mechanisms that modify the blockage effect in stably stratified flow.To investigate the impact from atmospheric stability, we simulate the flow around a wind plant for two distinct stable boundary layers.Furthermore, we run a set of simulations for a stand-alone turbine for each atmospheric condition to establish a baseline blockage effect for comparison.
This paper is structured as follows.We describe our simulation setup and stability cases in Sect. 2. We show how the velocity field in the induction region is modified by atmospheric static stability in Sect.3. In Sec. 4, we analyze the physical mechanisms modifying wind plant blockage.Finally, in Sect.5, we discuss our findings and provide suggestions for future work that could further improve understanding of the wind plant blockage effect.

Large-eddy simulation set-up
We perform LES of a wind plant under stably stratified conditions using the Weather Research and Forecasting (WRF) model v4.1.5(Skamarock et al., 2019) with turbines represented using a generalized actuator disk (GAD) approach (Mirocha et al., 2014;Aitken et al., 2014;Arthur et al., 2020).WRF is a fully compressible, nonhydrostatic model that solves the Navier-Stokes and thermodynamic equations for large-Reynolds number fluids (no viscosity or thermal conductivity).WRF uses an Arakawa-C grid staggering in the horizontal and a hydrostatic pressure-based vertical coordinate.Equations are integrated in time using a 3rd-order Runge-Kutta scheme with a smaller time step for acoustic modes.The advection terms are spatially discretized using a hybrid 5th-and 3rd-order scheme in the horizontal and vertical directions, respectively, which improves the model's effective resolution to 4 − 5∆x (Kosović et al., 2016).
We use a two-domain configuration with flat terrain to evaluate the blockage effect from wind plants.A periodic LES domain provides the boundary conditions for a nested LES domain via one-way nesting.Horizontal and vertical grid spacing remains the same for the parent and nest domains.We use the same domain characteristics, but a smaller nested domain size, to evaluate the blockage effect for a stand-alone wind turbine.We simulate moderately and weakly stably stratified flow.Because turbulence structures vary with atmospheric stability (Wurps et al., 2020), we use a finer grid in the more stable case.Horizontal grid spacing is ∆x = 5 m and ∆x = 7 m for the moderate and weak stability conditions, respectively.Similarly, vertical grid spacing at the surface is ∆z s = 5 m in the moderate stability case and ∆z s = 6 m in the weak stability case.A summary of the LES domains is in Table 1.All simulations are initialized with a dry atmosphere and zero latent heat flux.No cloud, radiation, or land surface models are used in the LES domains.An implicit Rayleigh damping term is applied to the vertical velocity in the upper 1000 m of 95 each domain to avoid wave reflection from the model top (Klemp et al., 2008).Surface boundary conditions are specified using Monin-Obukhov similarity theory with a surface roughness of z 0 = 0.1 m.We use the nonlinear backscatter and anisotropy (NBA) model with turbulence kinetic energy (TKE)-based stress terms (Kosović, 1997;Mirocha et al., 2010) to parameterize subgridscale (SGS) fluxes of momentum and heat.
Turbines are simulated exclusively in the nested domain using the generalized actuator disk implemented in WRF by Mirocha   (Jonkman et al., 2009).Both the wind plant and stand-alone turbine layouts used in our simulations are shown in Figure 1.Our wind plant has an aspect ratio of ∼ 3 : 2 to amplify the blockage effect as suggested by Allaerts and Meyers (2019).Segalini and Dahlberg (2020) found the blockage effect remains nearly constant when the wind plant has three or more rows, thus we include four turbine rows in our plant.Further, we constrain wind turbine spacing to 7D and 3.5D in the streamwise and cross-stream directions, respectively, following typical spacing in actual wind plants (Stevens et al., 2017).
Note that throughout the manuscript we denote temporal averaging using an overbar , spatial averaging along the i−direction using angled brackets ⟨ ⟩ i , and a normalized quantity using a hat ( ˆ).

Atmospheric conditions
We simulate two different boundary layers to evaluate how blockage varies with atmospheric stability (Figure 2).Distinct  A fully developed stable boundary layer is attained by spinning up a turbulent, neutral boundary layer, then adding a cooling rate at the surface.We initialize our simulations with a uniform potential temperature profile of θ = 300 K up to z = 1000 m, a capping inversion from 1000 m < z < 1200 m with ∂θ/∂z = 0.01 K m −1 , and we specify ∂θ/∂z = 0.001 K m −1 in the troposphere aloft.Both simulations are initialized with a U g = 11 m s −1 geostrophic wind speed.To reduce computational requirements, we spin-up turbulence and atmospheric stability in a small, precursor domain.After the flow is fully turbulent and stably stratified, we tile multiple precursor domains along the x− and y−directions to form a large domain.A complete description of this tiling methodology is presented in Appendix A.  A realistically turbulent neutral boundary layer develops shortly after initialization.Large, localized gradients of the horizontal velocity instigate large-scale turbulence early in the simulation, which then cascades into small-scale turbulence (Figure 3).Localized shear instabilities instigate turbulence throughout the boundary layer within the first hour of the simulation.These structures break up rapidly into smaller eddies, reducing shear until a quasi-steady state is reached.Turbulence structures form 125 rapidly close to the surface and propagate upwards (Figure 3).At hub height, turbulence propagates across all resolvable scales after 20 min (Figure 3a).Further aloft, large-and small-scale turbulent motions develop after 30 min (Figure 3b).Turbulence propagates up to the capping inversion after 50 min (Figure 3c).Even though the flow is fully turbulent 1 hr after initialization, turbulence and the mean flow in the boundary layer stabilize after 2 hr.Spin-up time varies for each stable simulation.The −0.5 K h −1 simulation is run until the bulk Richardson number in the surface layer and the Obukhov length become quasi-steady (Figure 6).The Obukhov length stabilizes to L = 35 m after 8 hr; the bulk Richardson number between z = 10 m and z = 153 m stabilizes to Ri bulk = 0.14.The nose of the LLJ after 8 hr is at z = 363 m for the −0.5 K h −1 case.The −0.2 K h −1 simulation is run for 13 hr.Even though L and Ri bulk do not reach a quasi-steady state (Figure 6), the flow displays characteristics typical of stable boundary layers and stability metrics suggest surface layer stability is different from the −0.5 K h −1 simulation.After 13 hr, the Obukhov length is L = 82 m and the bulk Richardson number is Ri bulk = 0.11.After 13 hr of simulation, the nose of the LLJ is at z = 660 m for the −0.2K h −1 case.Hub-height wind speed and direction after spin-up are comparable for both atmospheric conditions.For the moderate stability case, hub-height wind speed and direction are 8.9 m s −1 and 270.2 deg, respectively; for the weak stability case, 8.3 m s −1 and 269.9 deg, respectively.The winds at hub-height are in Region II of the power curve for the NREL 5MW turbine (Jonkman et al., 2009).
Each stability case is run for 1 hr, from which the first five minutes are discarded to guarantee that the flow at hub-height moves over the entire wind plant.The three-dimensional velocity, pressure, and potential temperature fields are output every 30 seconds, then time-averaged over the simulated time period.

The induction region
The three-dimensional, time-averaged velocity fields are averaged spatially to characterize the induction region.The velocity field is averaged vertically over the turbine rotor layer (27 m < z < 153 m).We also average the flow horizontally in the crossstream (i.e., y−direction) direction.We distinguish between the flow immediately upstream of each front-row turbine (i.e., intra-turbine) and the flow in between turbines (i.e., inter-turbine).Figure 7 illustrates the velocity deficit ∆U at hub height, defined as the difference between the time-averaged velocity at each grid cell U and the time-averaged velocity at the inflow of the domain U ∞ .The intra-turbine region is shown as the stippled area in Figure 7a,b.The inter-turbine region corresponds to the area in between the turbines (Figure 7a).The velocity deceleration immediately upstream (−2D < x < 0D) of the wind plant is strongly influenced by the induction region of the individual turbines (Figure 8).The velocity deceleration in the inter-and intra-turbine regions differs substantially within 2D upstream of the first row of the wind plant, but is roughly equal farther than 2D upstream (Figure 8).The flow decelerates more upstream of a wind plant compared to a stand-alone turbine (Figure 9).To compare the velocity deficit between the wind plant and stand-alone turbine directly, we consider the velocity in the intra-turbine region only.
On average for both atmospheric conditions, the velocity deceleration 2D upstream of the turbines is 1.57% and 3.27% for the stand-alone and wind plant, respectively.Likewise, the velocity deceleration extends farther upstream for a turbine array compared to a stand-alone turbine.Whereas the wind slows down by 1% relative to the domain inflow at 3D upstream of a stand-alone turbine, the same slowdown occurs 8.8D upstream of the wind plant.
Wind plant blockage is amplified with atmospheric stability (Figure 9).The flow upstream of the first turbine row decelerates more in the moderate stability case compared to the weak stability case (solid lines in Figure 9).Horizontal wind speed is on average 31% slower for the moderate stability case between x = −6D and x = −2D.For a stand-alone turbine, however, blockage does not change significantly for the stability regimes tested here (dashed lines in Figure 9).The velocity deficit is averaged vertically over the turbine rotor layer and horizontally over the inter-tubrine region (stippled area in Figure 7).The x−axis is scaled to locate x = 0D at the location of the turbine.
Even though the wind speed slowdown from blockage is small, front-row turbines in the wind plant produce on average 5.2% less power than a stand-alone turbine.Just as the velocity deceleration is modified with atmospheric stability, turbine underperformance is more severe in the moderate stability case compared to the weak stability case.Whereas turbines in the first row produce on average 4% less power than a stand-alone turbine for the weak stability condition, front-row turbines produce on average 6.5% less power than a stand-alone in the moderate stability case.

Physical mechanisms modifying blockage
Differences in upstream blockage for both stability cases can be explained by evaluating the forcing mechanisms driving the flow.The steady state integral momentum equation for the u−velocity (Eq. 1) provides insight into the physical mechanisms acting on the flow.In Eq. 1, vector quantities are in bold, î is the unit vector in the x−direction, and n is the outward pointing unit normal at each point on the surface S of the control volume V.
We evaluate the balance between momentum advection by the mean flow and a pressure gradient.Even though the Coriolis force in our simulation domain is not negligible, Coriolis forcing in the induction region is small.The Coriolis parameter scales as f c ∼ 10 −4 s −1 and the v−velocity in the turbine rotor layer is on the order of v ∼ 0.1 m s −1 , thus Coriolis forcing is of the order f c v ∼ 10 −5 m s −2 .Turbulence momentum redistribution is also small in the induction region of the wind plant In comparison, momentum advection by the mean flow in the induction region is of the order 10 −1 m s −2 .
We evaluate the integral momentum equation on differential control volumes to examine the streamwise evolution of the flow (Figure 10).Each differential control volume δV (blue rectangular cuboid in Figure 10) is bounded vertically within the  The balance between momentum advection and a pressure gradient along the x−direction for each differential control volume becomes: Each ∆(ρ u u i S i ) term represents the net advection of u−momentum by the u i velocity component through the S i control surface.For example, the net advection of u−momentum by the u−velocity in a control volume δV is ∆(ρ u u S x ) = (ρ u u S x ) out − (ρ u u S x ) in , as shown in Figure 10b.
The u−momentum flux at the inflow of the control volume V in Figure 10 is larger in the moderate stability case compared to the weak stability case due to slightly faster hub-height winds.Consequently, the magnitude of the momentum fluxes is expected to be larger in the moderate stability case as well.To contrast the momentum balance between different atmospheric stability conditions and turbine array sizes, we normalize the forcing terms in Eq. 2 using the momentum flux at the inflow of the control volume far upstream (ρ u ∞ u ∞ S x ) for each stability case.

Stand-alone turbine
Flow deceleration upstream of a stand-alone turbine is primarily caused by an adverse pressure gradient (Figure 11).The pressure gradient ∆pS x upstream of the turbine forms in response to the thrust force that the turbine imparts on the flow.The forcing mechanisms driving blockage for a stand-alone turbine remain virtually unchanged for the stability cases analyzed here (Figure 11).The pressure gradient force that drives flow deceleration upstream of the turbine differs by 3.1% between atmospheric conditions.Similarly, u−momentum advection by the v-and w−velocity components varies by 3.7% and 1.8%, respectively, between atmospheric conditions.

Wind plant
Momentum balance in the streamwise direction indicates flow deceleration upstream of the first turbine row is balanced by a pressure gradient and vertical momentum advection for both atmospheric conditions (Figure 12).An adverse pressure gradient ∆pS x forms immediately upstream of each front-row turbine, producing a force that decelerates the flow.The pressure gradient force decays rapidly upstream of the turbines.Also depleting momentum within the turbine rotor layer, the vertical advection of momentum ∆(ρ u w S z ) transports momentum upwards.Even though vertical momentum advection also decays upstream of the turbine array, the upwards transport of u−momentum remains larger than the pressure gradient force in the entire induction region of the wind plant for both atmospheric conditions.The cross-stream momentum advection ∆(ρ u v S y ) also reduces momentum availability, but only immediately upstream of the turbine array (−0.5D < x < 0D).The streamwise momentum advection ∆(ρ u u S x ) replenishes momentum upstream of the first turbine row.Vertical momentum advection is the primary forcing mechanism influencing wind plant blockage in our simulations.Figure 13 shows the net contribution of each term in Eq. 2 over the entire region upstream of the turbine array (control volume V in Figure 10a).Cumulatively over the induction region, vertical momentum advection is 41.3% (18.4%) larger than the pressure gradient force for the moderate (weak) stability case.Momentum advection by the v-velocity is 10.1% (12.8%) as large as the 240 vertical momentum advection for the moderate (weak) stability case.We now investigate how mean momentum transport and the adverse pressure gradient originate and compare within the induction region.

Pressure gradient force
Atmospheric stability marginally influences the streamwise pressure gradient upstream of the wind plant (Figure 14).The streamwise evolution of the normalized pressure gradient force remains nearly unchanged with atmospheric stability (Figure 245 14a).Furthermore, over the induction region, the normalized adverse pressure gradient upstream of the wind plant differs by 1% between the moderate and weak stability cases (Figure 14b).The pressure gradient upstream of a single front-row turbine in the wind plant is virtually the same as the pressure gradient force for a stand-alone turbine for both atmospheric conditions (Figure 15).The streamwise evolution of the pressure gradient does not vary significantly between −3D < x < 0D for a front-row turbine in the array and a stand-alone turbine (Figure 15a).Over the induction region of a stand-alone turbine (−6D < x < 0D), differences in the pressure gradient force between a turbine in the wind plant and a stand-alone turbine are smaller than 3% (Figure 15b).Given that the normalized pressure gradient force remains unchanged with atmoshperic stability and turbine array size, differences in blockage are caused by momentum redistribution in the induction region.

Mean momentum advection 255
The streamwise flow deceleration in the induction region is primarily transferred into upward motions (Figure 16).We evaluate  The vertical velocity advects horizontal momentum out of the turbine rotor layer (Figure 17).Vertical momentum advection is 20% larger in the moderate stability case compared to the weak stability upstream of the first turbine row (Figure 17b).
Larger vertical shear of the horizontal velocity in the moderate stability case contributes to the increased vertical momentum advection compared to the weak stability case.Shear (∆u/∆z) within the turbine rotor layer is 43.6% larger in the -0.5 K h −1 simulation compared to the -0.2K h −1 simulation.Similarly, the vertical velocity is 20% larger in the moderate stability case than in the weak stability case between x = −6D and x = 0D.Vertical momentum advection is amplified for a wind plant compared to a stand-alone turbine (Figure 18).For a given atmospheric condition, vertical shear of the horizontal velocity remains unchanged between the stand-alone turbine and wind plant simulations.Therefore, differences in vertical momentum transport are entirely due to the vertical velocity that forms upstream of the turbine array.For the wind plant, the secondary circulation (i.e., upwards w−velocity) extends farther upstream than for a stand-alone turbine.Whereas vertical momentum advection is on average 26% larger for a front-row turbine in the array compared to stand-alone turbine between −1D < x < 0D, it is x2 larger between −6D < x < −1D.Over the induction region of a stand-alone turbine (−6D < x < 0D), vertical transport of horizontal momentum is 72% (55%) larger for a frontrow turbine in the array than for a stand-alone turbine for the moderate (weak) stability case (Figure 18b).

Discussion and conclusions
The horizontal wind component within the rotor swept area decelerates upstream of wind plants, deflecting upward, a phenomenon called blockage.Blockage undermines energy production of wind plants by reducing turbine power production of mainly the front-row turbines.As the magnitude of the velocity deceleration and the nature of the physical mechanisms amplifying blockage are not yet well understood, we perform idealized WRF-LES of a finite-size wind plant and a stand-alone turbine for two atmospheric conditions to characterize blockage in stable boundary layers.Furthermore, we analyze the mechanisms driving and amplifying blockage by evaluating momentum conservation upstream of the turbines.
The velocity deceleration upstream of a wind plant is larger than upstream of a stand-alone turbine.For our simulations, the velocity deficit due to blockage is twice as large for the wind plant compared to a stand-alone turbine 2D upstream (Figure 9).Likewise, Reynolds-averaged Navier-Stokes simulations (RANS) by Bleeg et al. (2018) show the velocity slowdown 2D upstream of a wind plant is on average x1.9 larger than for a stand-alone turbine.For a variety of wind plant layouts, Strickland wind plant compared to a single turbine.For the simulations herein, the velocity deficit 7D upstream of a stand-alone turbine is negligible, whereas the deficit at the same distance for a wind plant is on average 1.2%.Comparably, Bleeg et al. (2018) suggests isolated turbines do not influence the flow 7−10D upstream, but wind plants do.However, it should be noted that Bleeg et al. (2018) do not provide a specific description of the stability cases or of the size of the wind plants used in their simulations, preventing a direct comparison.
The dominant physical mechanism that decelerates the flow upstream of a wind plant is different than for a stand-alone turbine in our simulations.Blockage for a stand-alone turbine is primarily caused by a pressure gradient upstream (Figure 11).The pressure gradient force is also present upstream of the wind plant; however, the primary mechanism decelerating the flow is the vertical advection of horizontal momentum (Figure 13).The pressure gradient force upstream of a front-row turbine of the wind plant is practically the same as for a stand-alone turbine (Figure 15).Conversely, the vertical transport of u−momentum upstream of a front-row turbine of the wind plant is on average 63% larger than for a stand-alone turbine between −6D < x < 0D (Figure 18).
Vertical advection of u−momentum is larger for a wind plant compared to a stand-alone turbine because of a secondary circulation that forms upstream of the turbine array.The slowdown of the u−velocity in the induction region of the wind plant is transferred into vertical motions (Figure 16).The vertical velocity advects horizontal momentum out of the turbine rotor layer.Given that the secondary circulation extends far upstream of the turbine array, vertical momentum advection is larger for the wind plant compared to the stand-alone turbine (Figure 18).
Boundary layer stability amplifies blockage for a wind plant (Figure 9).The wind speed is 3.5% and 2.8% slower than freestream 2D upstream of the wind plant for moderately and weakly stratified flow, respectively (solid lines in Figure 9).For a stand-alone turbine, atmospheric stability has only a second order effect on blockage.The wind slows down by 1.55% and 1.60% 2D upstream of a stand-alone turbine for the weak and moderate stability cases, respectively (dashed lines in Figure 9).Bleeg and Montavon (2022); Strickland et al. (2022) also found wind plant blockage changes with stability.Bleeg and Montavon (2022) quantify stability using the vertical change in potential temperature in the boundary layer ∆θ bl .For a 600 m deep boundary layer with weak temperature stratification (∆θ bl = 0.6 K), they report a single turbine row produces on average 2.2% less power than a turbine in isolation (Bleeg and Montavon, 2022).For a 600 m deep boundary layer with moderate stratification (∆θ bl = 2.9 K), they report a single turbine row produces on average 7.2% less power than a turbine in isolation (Bleeg and Montavon, 2022).In comparison, front-row turbines in our weak (z bl = 660 m, ∆θ bl = 1.37 K) and moderate (z bl = 360 m, ∆θ bl = 2.63 K) stability simulations produce 4% and 6.5% less power than a stand-alone turbine, respectively.
The differences in our results can arise from very different wind plant layouts, as blockage has been shown to be sensitive to turbine spacing and wind plant size (Strickland and Stevens, 2022;Bleeg and Montavon, 2022).Bleeg and Montavon (2022) perform RANS of 21 turbines in a single turbine row, while we do LES of a four-row wind plant with 10 turbines per row.
The fact that the velocity slowdown herein remains unchanged with atmospheric stability for a stand-alone turbine but not for The larger velocity deceleration upstream of a wind plant in the moderate stability case compared to the weak stability case is due to increased vertical momentum transport upstream of the first turbine row (Figure 13).The normalized pressure gradient in the induction region remains unchanged with atmospheric stability for the wind plant as a whole (Figure 14) and for individual turbines in the array (Figure 15).Conversely, vertical u−momentum advection upstream of the wind plant as a whole (Figure 17) and for individual turbines in the array (Figure 18) is larger in the moderate stability case compared to the weak stability case.Vertical shear of the horizontal velocity is 43% larger in the moderate stability case compared to the weak stability case, contributing to increased vertical momentum transport in the induction region.
Other studies analyzing the physical mechanism modifying blockage (in the absence of gravity waves) conclude the adverse pressure gradient amplifies blockage with closely spaced turbines (Strickland and Stevens, 2022;Strickland et al., 2022).Even though Strickland et al. (2022) demonstrate kinetic energy fluxes are the dominant mechanism modifying the power production of front-row turbines, they suggest an increased pressure gradient amplifies blockage as cold air is deflected upwards.
Note that we quantify the pressure gradient force upstream of the turbine array, which includes contributions from horizontal density gradients (i.e., cold air anomaly), and do not find significant changes with atmospheric stability (Figure 14).Strickland It is important to highlight that our simulations are idealized and in flat terrain.The role of terrain and vegetation should also be considered when evaluating wind plant blockage in future studies.Not only does terrain complicate efforts to measure blockage (Sanchez Gomez et al., 2022), but it likely also modifies the momentum balance upstream of the turbines (Segalini, 2017).Another important area of future research is whether or not wind plants trigger gravity waves, which may amplify   Even though the flow in our simulations is stably stratified, the effect of buoyancy is only significant in the moderate stability case (Figure B1).The vertical displacement of the flow in the induction regions is small, resulting in small changes in density over the induction region.For both stability cases, the streamline displacement at z = 153 m between the inflow of the domain and the first turbine row is smaller than 10 m.Consequently, the fractional change in fluid density in the induction region is 0.004% and 0.0002% for the moderate and weak stability cases, respectively.

Appendix C: Gravity waves
We examine the presence of gravity waves in our simulations by analyzing the correlation between the pressure, vertical velocity and potential temperature fields.We evaluate the streamwise deviation of each atmospheric variable from the inflow of we normalize each variable a i as âi = ai max(ai)−min(ai) , so that its values are between 0 and 1. Figure C1 shows the streamwise evolution of the deviation in vertical velocity, pressure and potential temperature from the inflow of the domain.There is no evidence of wind plant-triggered gravity waves in our simulation domain for either atmospheric condition.In our simulations, the pressure and vertical velocity are out of phase by ∼ 90 deg (Figure C1).Upstream of the wind plant, the vertical velocity perturbation reaches a maxima when the pressure perturbation is zero.At the outflow of the wind plant, the pressure perturbation is at a minima and the vertical velocity perturbation is close to zero.In internal gravity waves, the pressure perturbation is in phase with the vertical velocity perturbation (Banta et al., 1990).Furthermore, our simulations show the pressure perturbation and the potential temperature perturbation are out of phase by ∼ 45 deg (Figure C1).Conversely, the potential temperature perturbation and the pressure perturbation are 90 deg out of phase in gravity waves (Banta et al., 1990).

Figure 1 .
Figure 1.Relative location of the turbines in the wind plant (a) and stand-alone turbine (b) simulations for evaluating blockage.Forty NREL 5MW wind turbines constitute the 200MW wind plant simulated herein.
100 et al. (2014) and modified by Aitken et al. (2014) and Arthur et al. (2020).The NREL 5MW wind turbine has a hub-height of 90 m, a rotor diameter D of 126 m, cut-in speed at 3 m s −1 , rated speed at 11.4 m s −1 , and cut-out speed at 25 m s −1 stability conditions are obtained by providing different forcing at the surface.As suggested byBasu et al. (2008), we prescribe a cooling rate rather than a heat flux at the surface.Moderately and weakly stably stratified flow are obtained by forcing the boundary layer with Ṫs = −0.5 K h −1 and −0.2 K h −1 , respectively.

Figure 2 .
Figure 2. Wind speed (a,b) and potential temperature (c) profiles for the atmospheric conditions simulated herein.Atmospheric variables are averaged spatially over the entire domain and temporally over 1 hr.

Figure 3 .
Figure 3. Compensated turbulence spectra of the w−velocity for the ∆x = 7 m neutral LES at z = 90 m (a), z = 300 m (b), and z = 800 m (c).Spectra are color coded in 10-minute time increments since initialization.The dotted, black vertical line in each plot represents the effective grid resolution of our simulations.

Figure 4 .
Figure 4. Horizontal velocity profile for the ∆x = 7 m neutral LES.The velocity profile is averaged spatially over the entire domain.Profiles are color coded in 10-minute time increments since initialization.

Figure 5 .
Figure 5. Evolution of atmospheric variables for the moderate (a-e) and weak (f-j) stability simulations.Profiles are color-coded in 1-hr increments.

Figure 6 .
Figure 6.Temporal evolution of the Monin-Obukhov length (a) and bulk Richardson number (b) for each atmospheric state after prescribing a cooling rate at the surface.The bulk Richardson number is estimated between z = 10 m and z = 153 m.The shaded colored areas in each plot represent the simulation time for evaluating blockage.

Figure 7 .
Figure 7. Horizontal velocity deficit (∆U = U − U ∞) at hub height for the wind plant (a) and stand-alone turbine (b) simulations.The velocity fields are averaged in time over 1 hr of simulation.The stippled areas represent the intra-turbine regions.The inter-turbine regions are defined as the area in between turbines.

Figure 9 .
Figure 9. Normalized velocity deficit ∆Û = U −U ∞ U ∞ https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.turbine rotor layer, and horizontally in the y−direction by the area covered by the wind plant.Along the x−direction, each differential control volume is 15 m long.For the stand-alone turbine, the control volume is bounded in the y−direction by the rotor diameter (Figure 10c).

Figure 10 .
Figure 10.Illustration of the region considered in the analysis of the momentum balance along the x−direction for the wind plant (a) and stand-alone turbine (c).The integral momentum equation is evaluated on differential control volumes δV along the streamwise direction upstream of the wind plant (b).Each control volume is bounded vertically by the top (z = 153 m) and bottom (z = 27 m) of the turbine rotor layer.Horizontally in the y−direction, the control volume spans the region upstream of the wind plant (from y = 1953 m to y = 5922 m).For the stand-alone turbine, the control volume is bounded in the y−direction by the rotor diameter.Each differential control volume is 15 m long in the x−direction.The area of each control surface Si is shown in the differential control volume δV in Panel (b).
https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.Immediately upstream of the turbine (x = −0.1D), the pressure gradient force becomes negative because the GAD produces a pressure discontinuity on the flow.Momentum advection by the cross-stream ∆(ρ u v S y ) and vertical ∆(ρ u w S z ) velocity components also decrease momentum availability in the induction region of the turbine.Whereas the v−velocity transports momentum to both sides of the turbine, the w−velocity primarily transports momentum upwards.The streamwise velocity advects momentum back into the induction region of the turbine.

Figure 11 .
Figure 11.Streamwise evolution of the u−momentum equation (Eq.2) for a stand-alone turbine in the weak (a) and moderate (b) stability cases.The integral momentum equation is evaluated on differential control volumes δV along the x−direction, as shown in Figure 10.The x−axis is scaled to locate x = 0D at the location of the turbine.The mean momentum fluxes and the pressure gradient force are normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.

Figure 12 .
Figure 12.Streamwise evolution of the u−momentum equation (Eq.2) for the weak (a) and moderate (b) stability cases.The integral momentum equation is evaluated on differential control volumes δV along the x−direction, as shown in Figure 10.The x−axis is scaled to locate x = 0D at the location of the first turbine row.The mean momentum fluxes and the pressure gradient force are normalized using the momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.

Figure 13 .
Figure 13.Momentum balance over the entire induction region of the wind plant.The integral momentum equation is evaluated on the control volume V shown in Figure 10a.The control volume V is bounded in the x−direction by the inflow of the domain and the first turbine row (x = 5670 m).The mean momentum fluxes and the pressure gradient force are normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.

Figure 14 .
Figure 14.Streamwise evolution (a) and cumulative (b) pressure gradient force upstream of the wind plant.In panel (a), the integral momentum equation is evaluated on differential control volumes as shown in Figure 10a,b.In panel (b), the integral momentum equation is evaluated on the control volume V shown in Figure 10a.The pressure gradient force is normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.

Figure 15 .
Figure 15.Streamwise evolution (a) and cumulative (b) pressure gradient force upstream of a single turbine for each stability case.In panel (a), the integral momentum equation is evaluated on differential control volumes as shown in Figure 10c for a single turbine in the middle of the wind plant and for a stand-alone turbine.In panel (b), the integral momentum equation is evaluated on the control volume V shown in Figure 10c for a single turbine in the middle of the wind plant and for a stand-alone turbine.The pressure gradient force is normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.
the integral mass conservation equation ρ (u • n) dS = 0 on the differential control volumes shown in Figure10a,b.Mass balance indicates the slowdown of the u−velocity in the turbine rotor layer (∆(ρ u S x ) < 0) is balanced by the development of a secondary circulation in the form of net-upwards vertical motion (∆(ρ w S z ) > 0) for both stability conditions (i.e., ∆(ρ u S x )+ ∆(ρ w S z ) ≈ 0).The increase in vertical velocity is driven by a vertical pressure gradient, which balances the downward buoyancy force in the stably stratified flow (see Appendix B for a deeper analysis on vertical momentum balance).The change in v−velocity to the sides of the wind plant is only significant immediately in front of the turbine array (−0.5D < x < 0D).

Figure 16 .
Figure 16.Streamwise evolution of the mass conservation equation for the weak (a) and moderate (b) stability cases.The integral mass conservation equation is evaluated on differential control volumes δV along the x−direction, as shown in Figure 10.The x−axis is scaled to locate x = 0D at the location of the first turbine row.The mass fluxes are normalized using the mass flux at the inflow of the control volume far upstream (ρ u∞ Sx) for the respective stability case.

Figure 17 .
Figure 17.Streamwise evolution (a) and cumulative (b) vertical momentum advection upstream of the wind plant.In panel (a), the integral momentum equation is evaluated on differential control volumes as shown in Figure 10a,b.In panel (b), the integral momentum equation is evaluated on the control volume V shown in Figure 10a.The vertical momentum flux is normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.

Figure 18 .
Figure 18.Streamwise evolution (a) and cumulative (b) vertical momentum advection upstream of a single turbine for each stability case.In panel (a), the integral momentum equation is evaluated on differential control volumes as shown in Figure10cfor a single turbine in the middle of the wind plant and for a stand-alone turbine.In panel (b), the integral momentum equation is evaluated on the control volume V shown in Figure10cfor a single turbine in the middle of the wind plant and for a stand-alone turbine.The vertical momentum flux is normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx) for the respective stability case.
https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.and Stevens (2022); Strickland et al. (2022) also show the velocity deceleration 2D upstream is consistently larger for a turbine array compared to a single turbine in isolation.The velocity deficit in the induction region also extends farther upstream of a https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.a wind plant suggests differences in blockage for the wind plant are due to the large-scale interaction between the turbine array and the boundary layer.

and
Stevens (2022);Strickland et al. (2022)  do not consider the influence of momentum advection upstream of the turbines; nonetheless, they show increasingly larger vertical velocity upstream of the turbines with smaller cross-stream turbine spacing and stronger atmospheric stability, likely increasing momentum advection outward of the turbine rotor layer.Furthermore,Strickland and Stevens (2022)  simulate purely neutral flow whereas we simulate stably stratified flow in a boundary layer with a capping inversion above.Strickland et al. (2022)  simulate an infinitely wide wind plant without a capping inversion.Bleeg and Montavon (2022) show the velocity reductions from blockage, and thus the amplifying mechanisms, vary when a capping inversion, a typical feature of the planetary boundary layer, is not included in the simulations.Thus, while wind plant blockage is caused by an adverse pressure gradient forming upstream the turbines, it is amplified by the vertical advection of momentum in the turbine rotor layer for the stability cases analyzed here.A pressure gradient forms from the thrust force of the individual turbines.The cumulative deceleration of the flow upstream of the individual turbines is balanced by an increase in vertical velocity.Vertical motions transport horizontal momentum out of the turbine rotor layer, reducing momentum availability at the turbine locations.
blockage.Future field experiments should evaluate gravity-wave initiation around wind plants.Knowledge about interactions between gravity waves and wind plants may provide a better understanding of the physics of wind plant blockage and aid in model validation.https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.

Figure A2 .
Figure A2.Correlation of atmospheric variables between adjacent tiles in the large LES domain for the weak (top) and moderate (bottom) stability cases.Profiles are color-coded in 5-minute increments since initialization of the large domain.The grey, shaded region represents the turbine rotor layer.

Figure B1 .
Figure B1.Streamwise evolution of the w−momentum equation (Eq.B1) for the weak (a) and moderate (b) stability cases.The integral momentum equation is evaluated on differential control volumes δV along the x−direction, as shown in Figure 10.The x−axis is scaled to locate x = 0D at the location of the first turbine row.The mean momentum fluxes, the pressure gradient force and buoyancy are normalized using the u−momentum flux at the inflow of the control volume far upstream (ρ u∞u∞ Sx).

Figure C1 .
Figure C1.Streamwise evolution of the vertical velocity, pressure and potential temperature deviation from inflow conditions for the weak (a) and moderate (b) stability cases.Each variable ai is normalized as âi = https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.or the U.S. Government.The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow 435 others to do so, for U.S. Government purposes.JDM and RSA contributed under the auspices of the U.S. Department of Energy (DOE) by Lawrence Livermore National Laboratory, under contract DE-AC52-07NA27344.https://doi.org/10.5194/wes-2023-20Preprint.Discussion started: 21 February 2023 c Author(s) 2023.CC BY 4.0 License.