A Subgrid Parameterization for Wind Turbines in Weather Prediction Models with an Application to Wind Resource Limits

A subgrid parameterization is offered for representing wind turbines in weather prediction models. The parameterization models the drag and mixing the turbines cause in the atmosphere, as well as the electrical power production the wind causes in the wind turbines.The documentation of the parameterization is complete; it does not require knowledge of proprietary data of wind turbine characteristics. The parameterization is applied to a study of wind resource limits in a hypothetical giant wind farm.The simulated production density was found not to exceed 1Wm, peaking at a deployed capacity density of 5Wm and decreasing slightly as capacity density increased to 20Wm.


Introduction
Wind power production in numerical weather prediction models can be either inert or active.In the inert type, the wind speed forecasted for a turbine location can be extracted from the model and used to calculate wind power production, with no impact of the turbines on the weather prediction [1].In the active type, this impact is included, specifically the drag and turbulence enhancement of the wind turbine acting on the atmosphere [2].In this paper we offer some details of a wind turbine parameterization appropriate for large wind farms, with many turbines within a grid cell.This paper refines the wind turbine parameterization in [2,3], effectively offering a simplified and documented alternative to what appeared in WRFv3.3 [4,5].Being subgrid, wakes are not explicitly simulated, but rather the momentum loss is immediately diffused across the breadth of the grid cell.The parameterization is adaptable to typical wind turbine characteristics.The giant wind farm of [3] is revisited for the purpose of studying the practical limit to wind power extraction from the atmosphere.Whereas [3] examined the much more subtle effect of the wind farm on precipitation climate statistics, the current study is more straightforward and does not require multidecadal simulations.The simulations use WRFv3.1 with the MYJ boundary layer scheme and 30 km horizontal grid spacing.The wind turbine parameterization adds elevated drag and production of turbulent kinetic energy to the MYJ scheme.
If a horizontal wind vector ⃗  is known at the height of wind turbine (in practice, meaning that a suitable average wind vector is known), then the power produced is where   is the capacity coefficient and  max is the rated power output for the particular wind turbine.In the simulations, we focus on  max = 2 MW and  max = 8 MW, which roughly brackets the range of potential installations.  is constrained by both laws of nature and engineering design.For  <  in , the turbine blades do not rotate, so   = 0. Likewise, for  >  out the turbine rotation is halted to avoid damage, and   = 0.As  increases past  in , power production rises rapidly, but by engineering design is brought to a broad plateau of  max , by mechanical adjustment of the turbine blade pitch angle [6].
Figure 1 shows a typical   () offered by this parameterization.Figure 1  on the atmosphere is to be calculated.The aerodynamical basis of (1) determines those impacts.From elementary physics, the available power of wind impinging on the rotor cross-sectional area  of the turbine is where  is the density of the air.The ratio of  to the available power is the power coefficient   : So (1) can be written as Though (4) provides the identical calculation of power production as (1), knowledge of   will have another important purpose in calculating the production of turbulent kinetic energy.
The drag force ⃗  on an object presenting cross-sectional area  to a uniform stream of fluid with velocity ⃗  is conventionally modeled in terms of a shape-dependent drag coefficient   .In the language of wind turbine modeling, the drag coefficient is named the thrust coefficient   : At large Reynolds number,   is predominantly shape dependent.For example, there are many references giving values such as for a flat plate   = 1.28 and for a sphere   = 0.47.The drag force for rotating turbine blades is much greater than a calculation based on stationary blades and using just the area presented by the blades.The drag force of the wind turbine is characterized in terms of the disk swept out,  =  2 , where  is blade length.For the Bonus Energy A/S 2.0 MW,   peaks at approximately   = 0.88 [2].Presumably this cited value of   includes the drag of the tower as well, contained between two pressure levels (green).In this example, the hub of the rotor lies between the pressure levels.The area above the hub is sum of two triangles and two sectors.The area of the triangles sum to ℎ √  2 − ℎ 2 .The area of the two sectors sum to  2 sin −1 (ℎ/) In this example, the area below the hub makes a similar positive contribution.If the hub is not between the layers, the total area is given by the subtraction of two areas calculated from the hub.Likewise, if ℎ > , then ℎ is replaced by  in the calculation.
but in this parameterization the drag force is modeled as occurring within the area of the rotor.
The drag force of the wind turbine removes momentum from the atmosphere and transfers it to the Earth.But with the Earth having a large mass, the drag force transmitted to Earth, via the tower, does not do significant work on the Earth, meaning that the loss of energy from the mean wind goes into power production and turbulent kinetic energy, rather than into kinetic energy of the Earth [4].
The force of the turbine on the atmosphere is opposite to that of (5), so the rate of work (power)   on the atmosphere is − ⃗  ⋅ ⃗ : By our principle of strict energy conservation so Most numerical weather prediction models employ a force per mass at a grid point, within a grid volume.The drag force in (5) would need to be normalized appropriately, by the total mass of air in the grid volume.Similarly, a normalization is required when ( 8) is used to predict turbulent kinetic energy and added to the other source terms in the prediction for turbulent kinetic energy.
Here we take  as the only portion of the multiple wind turbine areas that are within the heights bounding a grid volume (Figure 2).This introduces some unrealism, as it allows a wind turbine to be modeled as having different rotation speeds and different   at various heights.The normalization procedure of the previous paragraph means that the rotor area per grid volume (an area density with units of inverse length) is the quantity needed for computation.In a staggered grid model, the heights of the prediction of horizontal wind may lie between the levels for the prediction of vertical velocity, the levels of which define the vertical bounds to the grid volume for horizontal velocity.

Functions for 𝐶 𝑓 (𝑉) and 𝐶 𝑡 (𝑉)
For   (), we employ a soft-clip function, which allows   () to come to a plateau without a sharp "knee point." Note that [4] does not provide this soft-clip feature and [2,3] do not have a monotonic   .The soft-clip function that we employ is computationally efficient and provides a very close approximation to (1 + tanh())/2: Here  = ( −  0 ), where  0 controls the center point of   and  controls the slope of the transition.Adjusting the center point and slope to approximate the characteristics of the Bonus wind turbine is elementary.Slightly more complicated is to require   to be exactly zero for  <  in , which may require employing a shift of () to bury part of it below the -axis.Let Thus As the blades of a wind turbine are adjusted, to reduce   so that the maximum in   does not exceed 1, the thrust coefficient is also reduced.We find the following fit satisfactory: Figure 1 lists the values of the parameters used for the Bonus 2 MW turbine.We employ a value for  that makes   ( out ) =   , so there is no discontinuity in   at  out .
Other choices are possible.

Application to Wind Resource Limits
Textbooks in atmospheric science cite the typical midlatitude pressure gradient force (per mass) to be 10 −3 m s −2 and the typical horizontal velocity scale to be 10 m s −1 , with about  Theses plots show the details behind the points for CD = 2.5 W m −2 in Figures 7 and 8.The average production with 8 MW turbines is 106 GW.The average production with 2 MW turbines is 66 GW.
1/10 of that being cross-isobaric.Only the cross-isobaric component is capable of internally renewing kinetic energy that has been removed by the wind farm.In the boundary layer, a larger fraction of the wind vector could be crossisobaric, but the magnitude of the vector could be less.So if we accept 1 m s −1 as the typical magnitude of crossisobaric flow, the rate of kinetic energy production within the wind farm would be 1 W m −2 per kilometer of depth of the extraction (assuming a density of  = 1 kg m −3 ).
Kinetic energy also enters the wind farm from the side.If a giant wind farm of horizontal area  ×  is extracting wind from a layer of depth  (where  could be the depth of the atmospheric boundary layer, rather than height to the top of the wind turbine), then (2) can be used to calculate the power advected into the wind farm presenting a side area of  =  × .This is power that can be potentially extracted over the area of the wind farm, giving a power density

Advances in Meteorology
For example, let us take  = 1 km and  = 100 km.For / = 0.01 and  = 5m s −1 , that gives an upper limit to extraction of 0.625 W m −2 .Using  = 10 m s −1 instead gives 5 W m −2 .Note that a giant wind farm could have / < 0.01, and the bound on power that is extractable from the advection source would correspondingly be less.Thus as / becomes small, renewing of the wind resource by the pressure gradient becomes more important.
In the above estimates, what value should be used for ?Also, what is the contribution of transport through the top of the wind farm?We also need to recognize that as the power extraction approaches the upper limit, that would imply that  would be decreasing as the wind farm is traversed.The continuity equation would thus require upward advection of energy out of the top of the wind farm.All these considerations imply that a more refined estimate of the limits to power extraction at a site will require details about the wind climate, including boundary layer mixing, as well as the use of a numerical weather prediction model.
Here we demonstrate an application of our wind farm parameterization with a modest extension to several such studies of limits to wind farm resources [2,4,7], namely, allowing for all sizes of wind turbines to have the characteristics of Figure 1.The parameterization is used to investigate the relative effect of deploying the capacity density as 8 MW turbines, twice the size of 2 MW turbines.
The wind farm location is as in [3], with an area of 182,700 km 2 (Figure 3).In [3], 2 MW wind turbines were situated with turbine density of 1.25 km −2 , giving the giant wind farm a capacity of 457 GW and a capacity density of 2.5 W m −2 .The 2 MW wind turbines had a hub height of 60 m and a rotor radius of 38 m.Here we experiment with both 2 MW and 8 MW turbines, deployed with capacity density ranging from 0.625 W m −2 to 20 W m −2 .The 8 MW wind turbines are simply double in height and radius of the 2 MW turbines.The   ,   , and   curves for both models are as in Figure 1.The study shown here is much simpler than [3], and examines only the effect of the wind farm on wind characteristics within the wind farm, as well as the power production by the wind farm.
The 8 days from 0 UTC April 23, 1948 to 0 UTC May 1, 1948 were convenient for this study.The National Renewable Energy Laboratory displays the annual average wind speed at 80 m to range between 7.0 m s −1 and 9.0 m s −1 in the modeled wind farm area.The model, without wind turbines, has an average wind speed of 7.0 m s −1 and 8.1 m s −1 at 52 m and 102 m, respectively, during the 8 days of the simulation.

8 MW versus 2 MW Deployment.
Here we highlight a particular comparison between deploying the 457 GW as either 228,375 2 MW turbines or as 58,656 8 MW turbines.In the analysis of power production, the average   experienced over the entire farm is denoted by CF (  is an engineering design parameter and CF is an experimental result).Though the capacity density (CD) is the same, the production density (PD) with the 8 MW deployment is 60% greater (Figure 4).A naive estimate might have anticipated an increase greater than 100%, using reasoning that the layer being mined for power is twice as deep, with the upper part having stronger winds.That sort of estimate is not realized.
Figures 5 and 6 show that the extraction (as indicated by wind speed reduction) has become rather insignificant at height 648 m.But both the 8 MW and 2 MW turbines discernibly remove energy below 648 m.As expected, the taller 8 MW turbines extract more energy in the layer above the height of the smaller turbine.Estimating how much more effective the extraction is with taller turbines has required the benefit of a numerical simulation.The ability to make such an estimation is one of the main practical benefits of the parameterization.

Production Saturation.
Here we summarize the investigation into power production across a broad range of wind farm characteristics.Conclusions are similar to [2,4,7,8]: Figure 7 shows a limit to power extraction to be on the order of 1 W m −2 .Such knowledge obviously influences design characteristics of wind farms: whether to add more wind turbines to a farm, acquire more land and develop a larger farm, or develop another farm in a distant location.In our simulations, the drop in CF proceeds immediately from the lowest CD.This is because power production is very sensitive to changes in  in the vicinity of  in , with   rising faster  than  3 wherever   is increasing with .The opposite scenario could happen in a different wind climate.If  is consistently well into the range that produces   () = 1, there might be a significant decrease in , but no drop in CF until CD exceeds a value greater than 1 W m −2 .
Consider increasing CD from 2.5 W m −2 , with 2 MW turbines, to 10 W m −2 by either quadrupling the area of the rotors or quadrupling the number of turbines.The two scenarios can be found within Figure 8. Increasing the rotor area density by a factor of 4 (redeploying as 8 MW turbines) increased PD by 2.07.Quadrupling the number of 2 MW turbines increased PD by a factor of 1.38.Since the increase in PD was significantly less than 4, we would say that the collective impact of the turbines on the power productivity of the winds is significant.Inspection of the wind difference plots in Figure 6 shows wind being reduced above the tops of the turbines, evidently the effect of turbulent transport of momentum vertically in the atmosphere.This transport may be hard to estimate by means other than a detailed numerical model.
We note that Horns Rev 1, an established 20 km 2 wind farm in the North Sea, has been averaging PD = 3.98 W m −2 in the last 5 years [9].This illustrates the importance of scale in understanding the production limitations of wind energy.As discussed in Section 3, the larger the horizontal extent of the wind farm is, the less important the advection of kinetic energy is and the more important the pressure gradient force becomes in sustaining energy producing wind speeds within the wind farm.The Horns Rev 1 wind farm is small enough that / is approximately 0.2, thus explaining the observed PD.

Conclusions
When considering national and international energy portfolios, wind energy continues to become an important part of diversified energy portfolios.Though current wind farms are small enough in scale to have / ratios that allow advection of kinetic energy into the side of a wind farm to be an important power source, it is important to discern how that wind resource diminishes with larger wind farms.Future power needs could force the development of giant wind farms, with areas that are orders of magnitude larger than current farms.Furthermore, the development of many small wind farms in close proximity could have a resource limit similar to a giant wind farm.
Giant wind farms will need to be planned with an active type numerical weather prediction model, so as to get an accurate estimate of the wind power resource.For example, in our study, larger (taller) wind turbines produce a larger CF.The cost-effectiveness of deploying larger turbines will require an accurate prediction of this CF before a financial decision can be made.

Figure 2 :
Figure2: Calculating the fraction of the area of the rotor circle (red) contained between two pressure levels (green).In this example, the hub of the rotor lies between the pressure levels.The area above the hub is sum of two triangles and two sectors.The area of the triangles sum to ℎ √  2 − ℎ 2 .The area of the two sectors sum to  2 sin −1 (ℎ/) In this example, the area below the hub makes a similar positive contribution.If the hub is not between the layers, the total area is given by the subtraction of two areas calculated from the hub.Likewise, if ℎ > , then ℎ is replaced by  in the calculation.

Figure 3 :
Figure 3: The wind farm is within the green rectangle.The average wind difference at 102 m between the simulation with 2.5 W m −2 of 2 MW turbines minus the simulation without turbines is shown.

Figure 4 :
Figure 4: The average wind farm capacity factor CF for two particular deployments of CD = 2.5 W m −2 , as a function of time.Theses plots show the details behind the points for CD = 2.5 W m −2 in Figures7 and 8.The average production with 8 MW turbines is 106 GW.The average production with 2 MW turbines is 66 GW.

− 1 )Figure 5 :
Figure 5: Average horizontal wind speed at the lowest 8 grid wind levels of 15, 52, 102, 165, 245, 345, 466, and 648 meters above ground.Orange is without wind turbines.Blue is with 2 MW wind turbines deployed at CD = 2.5 W m −2 capacity.Top of the rotor is at 98 m (careful counting shows that blue curve is below the corresponding orange curve).The wind barbs represent the horizontal wind direction and magnitude at 102 m at the center of the wind farm area, but without wind turbines.Half barbs are 1 m s −1 , full barbs 2 m s −1 , and flags 10 m s −1 .Days are ticked at 0 UTC, late afternoon at the wind farm, at which time the wind speed has become more well mixed across the boundary layer.

Figure 6 :Figure 7 :
Figure 6: Area-averaged wind speed change that results from deploying wind turbines, as a function of height and time.The height of the top and bottom of the rotor are indicated by longer tick marks at the extreme right.The extent of the wind reduction above the tops of the turbines is indicative of power extraction from those layers, a complicated prediction requiring a numerical weather prediction model.

2 )Figure 8 :
Figure8: Average wind farm production density PD for various deployments.The two points marked with a star are for simulations repeated using twice the number of grid points within the depth of the boundary layer, as compared with the standard resolution.The simulations with standard resolution are indicated with a circle.The vertical positions of the grid points in the standard resolution are listed in Figure5.