From shear to veer: theory, statistics, and practical application

. In the past several years, wind veer — sometimes called ‘directional shear’ — has begun to attract attention due to its effects on wind turbines and their production, particularly as the length of manufactured turbine blades has increased. Meanwhile, applicable meteorological theory has not progressed significantly beyond idealized cases for decades, though veer’s effect on the wind speed profile has been recently revisited. On the other hand the shear exponent ( α ) is commonly used in wind energy 5 for vertical extrapolation of mean wind speeds, as well as being a key parameter for wind turbine loads calculations and design standards. In this work we connect the oft-used shear exponent with veer, both theoretically and for practical use. We derive relations for wind veer from the equations of motion, finding the veer to be composed of separate contributions from shear and vertical gradients of cross-wind stress. Following from the theoretical derivations, which are neither limited to the surface-layer nor 10 constrained by assumptions about mixing length or turbulent diffusivities, we establish simplified relations between the wind veer and shear exponent for practical use in wind energy. We also elucidate the source of commonly-observed stress-shear misalignment and its contribution to veer, noting that our new forms allow for such misalignment. The connection between shear and veer is further explored through analysis of one-dimensional (single-column) Reynolds-averaged Navier-Stokes solutions, where we confirm our theoretical derivations as well as the dependence of mean shear and veer on surface roughness 15 and atmospheric boundary layer depth in terms of respective Rossby numbers. Finally we investigate the observed behavior of shear and veer across different sites and flow regimes (including forested, offshore, and hilly terrain cases) over heights corresponding to multi-megawatt wind turbine rotors, also considering the effects of atmospheric stability. From this we find empirical forms for the probability distribution of veer during high-veer (stable) conditions, and for the variability of veer conditioned on wind speed. Analyzing observed joint probability distributions of 20 α and veer, we compare the two simplified forms we derived earlier and adapt them to ultimately arrive at more universally applicable equations to predict the mean veer in terms of observed (i


Introduction
The shear exponent has generally not been used or accepted by meteorologists, as it does not (directly) relate to the physics of atmospheric flow, nor to the most important boundary condition-the surface.Regarding the latter, in contrast with similarity theory (Monin and Obukhov, 1954), the shear exponent does not contain explicit information about the surface roughness.
However, the shear exponent can be related to surface properties in a generalized way, as well as to turbulent kinetic energy and atmospheric stability (buoyancy) as shown by e.g. Kelly et al. (2014a).This is particularly useful above the atmospheric surface layer (ASL), where micrometeorological theory based on ASL assumptions fails-and where the effects of the surface are neither dominant nor simple enough to be characterized through accepted ASL parameterizations.As practiced in the wind energy resource assessment community for decades, the shear exponent can thus be preferable over similarity theory for use in vertical extrapolation (Irwin, 1979;Mikhail, 1985;Petersen et al., 1998) with quantification of uncertainty in its use more recently reinforcing such (Triviño et al., 2017;Kelly et al., 2019b).Shear is also a key parameter for flow characterization towards loads simulations, being seen to systematically affect various turbine loads (e.g.Dimitrov et al., 2018;Robertson et al., 2019).
Veer has received much less attention than shear, though its potential importance to wind energy has been noted more recently.In the meteorological literature, where veer is often labelled as 'directional shear' or 'turning, ' Markowski and Richardson (2006) reviewed the distinction between veer and vertical gradients of wind speed, listing studies of meteorological phenomena that considered veer (though they focused on convective storms).While some works in meteorology have investigated veer, these have tended to focus on the angular difference between winds at the top of the atmospheric boundary layer (ABL) and the surface (e.g.Clarke, 1975;Brown et al., 2005;Grisogono, 2011;Lindvall and Svensson, 2019),and are not generally suited for engineering applications.For wind energy, Murphy et al. (2020) looked at the veer (and shear) along with power production measured over a six-month period, finding a minor but non-negligible effect of veer on power production for a utility scale turbine.Gao et al. (2021) found positive veer over the upper half of a single (2.5MW) clockwise-turning turbine rotor to reduce power production, opposite and slightly larger than the corresponding effects of negative veer there; they also showed the rotor's lower-half veer was less significant than the upper half.Shu et al. (2020) examined measurements from a lidar offshore between islands southwest of Hong Kong, observing larger veer when hilly terrain was upstream compared to more open sea conditions; they also noted seasonal variations.For power production, the veer was incorporated into rotor-equivalent wind speed (REWS) by Choukulkar et al. (2016), whom found it to generally decrease production at two sites; Clack et al. (2016) found similar results from weather assimilation model output over the USA, along with higher production at night and lower power during daytime at most locations.Wind veer has also been examined with regard to its connection with the distortion and lateral movement turbine wakes via measurements and simulations (e.g.Abkar et al., 2018;Brugger et al., 2019), also including yaw-misalignment affects (Hulsman et al., 2022;Narasimhan et al., 2022).
In this paper we elucidate analytical and statistical connections between a number of key parameters used to describe atmospheric boundary layer flow, with focus on the vertical variation of wind velocity.This follows the earlier work of Kelly et al. (2014a) that gave forms for low-order statistics of shear exponent α, relating α to turbulence intensity and stability; here we derive new relations for the turning of the wind in terms of shear, going beyond classical Ekman-type analysis.
2 Theory and development

Shear exponent
Just as potential temperature-the buoyancy variable commonly-used in meteorology-was labeled the "meteorologist's entropy" by Bohren and Albrecht (1998), one could call the shear exponent (α) the "wind engineer's phi-function."Specifically this follows from the definition of shear exponent and the dimensionless wind speed gradient used in meteorology, where u * 0 is the surface layer friction velocity (square root of kinematic shear stress), κ = 0.4 is the von Kármán constant, and z is the height coordinate 1 .We remind that (1) is derived from the power-law expression for wind speed which is assumed to be valid over some extent around height z ref , with U ref ≡ U (z ref ).The power-law (3) with shear exponent (1) has been used in wind engineering for decades (e.g.Irwin, 1979;Petersen et al., 1998) due to its simplicity, and because it doesn't require any information other than the wind speed at two heights.Although (1) and ( 2) might appear to be quite alike, one can see a phenomenological difference when comparing the wind speed profiles resulting from these relations.In Monin- is the M-O wind speed correction function; the analytic forms for Φ m and Ψ m differ in stable and unstable conditions, and have been determined empirically in decades past (Businger et al., 1971;Carl et al., 1973;Li, 2021).But Monin-Obukhov similarity theory and its assumptions (such as constant u * ), as well as established forms for Φ m , fail above the surface layer;3 this motivates use of α in applications such as wind energy, as (1) does not directly rely on surface-layer assumptions.

Relation to stability and turbulence
As shown by Kelly et al. (2014a), in horizontally homogeneous conditions the steady or mean balance of turbulent kinetic energy (TKE) can be written in terms of shear exponent as for a given height z, where the streamwise direction is defined by the mean wind U (z) and we have suppressed z-dependences for brevity; here ⟨uw⟩ is the turbulent horizontal momentum flux (kinematic stress), T is the total (turbulent plus pressure) transport, B is buoyant production, and ε is the viscous dissipation rate of TKE.Within the ASL under these conditions where M-O theory is valid, using the neutral value of dissipation rate as ε 0 ≡ u 3 * 0 /(κz) along with the dimensionless functions Φ ε ≡ ε/ε 0 and Φ T ≡ T /ε 0 (Kaimal and Finnigan, 1994), we can express an ASL version of (4) as since by definition B = −z/L and u 2 * 0 = −⟨uw⟩; here I u ≡ σ u /U is the streamwise turbulence intensity.The dimensionless dissipation rate (M-O function) Φ ε ≥ 1 is roughly 1 + 5z/L in stable conditions and increases more weakly with −z/L in unstable conditions (Panofsky and Dutton, 1984;Kaimal and Finnigan, 1994); meanwhile the transport is negligible in stable conditions but Φ T > 0 in unstable conditions (e.g.Wyngaard, 2010).Thus in stable conditions (L −1 > 0) one can see α is larger than in neutral conditions, while in unstable conditions α becomes smaller.Above the ASL this will also generally be the case, though analytic nondimensional forms become difficult to derive, while the flow becomes affected by more terrain upwind and associated inhomogeneities; furthermore in stable conditions the local stability (at a given z) becomes increasingly more important than surface-based z/L (Derbyshire, 1990).As will be shown below, the most common and mean conditions at contemporary rotor heights qualitatively follow ( 5), but due to these and other non-ideal effects (e.g.nonstationary transients) large deviations can occur.We note that in this work we are not searching for analytical forms for α or surface-layer behavior; rather, we are concerned with how α relates to the veer, especially over heights corresponding to wind turbine rotors, a portion of which commonly extends beyond the ASL.https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.

Veer
For the simplified general case of Coriolis-affected mean flow, we write the horizontal mean velocity vector {U, V } as a complex number, S ≡ U + iV = |S|e iφ .For a mean wind direction defined at some height z, the veer can be defined as a directional shear ∂φ/∂z through the wind direction In most of the micrometeorological literature, the mean wind direction is defined based on the surface stress (i.e. via the winds closest to the surface, so φ 0 ≡ φ(0) = 0).We follow this convention unless stated otherwise, as done for some expressions later in section 2.3; one could also choose to define the coordinate system based on the geostrophic wind direction (e.g.Svensson and Holtslag, 2009).
As is classically known in micrometeorology (e.g.Hess and Garratt, 2002), the veer across the entire ABL depends primarily on the Coriolis parameter f (thus latitude), geostrophic wind speed |G|, and surface roughness length z 0 , but is also affected by the ABL depth h and stability (as confirmed via Reynolds-averaged Navier-Stokes simulations by van der Laan et al., 2020).
The veer across a fraction ∆z/h of the ABL will also depend on these parameters; thus for a given site and height, ∆φ/∆z will have a distribution due to variations in these parameters.This will become clearer below as we examine the relationship between veer and shear.
The Coriolis-affected mean momentum balance can be written in the form for stationary and horizontally homogeneous conditions (thus neglecting advection).Here the kinematic horizontal pressure The mean stresses are dominated by vertical momentum transport ⟨sw⟩, where w denotes (turbulent) vertical velocity fluctuations and s ≡ u + iv the horizontal velocity fluctuations.
At a given height z, taking the differential of (6) (recalling d arctan x = dx/[1 + x 2 ] and using the chain rule) gives here the superscript asterisk denotes complex conjugate.Applying ∂/∂z to ( 8) and ( 7) and combining provides a basic expression for veer: In the case of zero geostrophic shear (dG/dz = 0), if the coordinate system's x-axis is defined by the mean wind direction at the height z where the veer is sought, then (9) can be written more simply as Though ( 9) and ( 10) are not directly very useful for relating veer to shear, they illustrate that the curvature of stress profiles and Coriolis effect are the basis for mean veer following (7), and also that geostrophic shear can further contribute to veer (e.g. due to baroclinity, Hoxit, 1974;Arya and Wyngaard, 1975;Pedersen et al., 2013).

Relating veer to shear
Towards relating the veer to shear one can alternately take the time derivative of (8); using the real and imaginary parts of ( 7), in the horizontally homogeneous limit (ignoring advection) one obtains a rate equation for mean wind direction: The 'turning' angle γ ≡ φ − φ G between geostrophic and mean wind directions (e.g.Wyngaard, 2010) arises through 4 by taking ∂/∂t of (6) or equivalently U ∂V /∂t−V ∂U/∂t via (8).The geostrophic wind direction is defined as and the 'cross-isobar' angle, i.e. the turning over the whole ABL (γ 0 = φ 0 −φ G ), is generally less than 45 • (Grisogono, 2011) 5 ; in a right-handed coordinate system, regardless of whether x is chosen to align with G or the surface-layer wind velocity U ASL , the turning tends to γ > 0 in the Northern hemisphere 6 .We remind that φ, and thus γ, can vary with height z (as can φ G in baroclinic conditions).
Assuming statistical stationarity so that ∂φ/∂t = 0, the vertical derivative of (11) can be written most conveniently in terms of the dimensionless deviation of the wind from streamwise; taking the vertical derivative of (11) if we again take dG/dz = 0 (neglect baroclinity), then As it is expressed in terms of angular differences γ, the equation above is independent of whether the coordinate system is defined at the surface or by the geostrophic wind.Expression (12) clearly separates the shear and Coriolis/stress contributions to veer.However, it can be simplified, and is most meaningful, if the coordinate system is defined at the height z for which it is applied; in practice the veer is typically calculated around hub height, or from hub to tip, or between measurement and hub heights.Re-expressing (11) with the coordinate system defined by having x in the mean wind direction at height z, so that 4 The turning angle can also be expressed in complex notation, recalling that the angle between vectors written in complex notation (here U → S and G → G) can be recovered by taking ℜ{G * S}, i.e. |G||S|ℜ{e −i(φ+γ) e iφ } = |G||S| cos γ. 5 The ABL turning angle γ 0 cannot exceed 45 • , according to the Ekman equations (or their numerical solution, as in van der Laan et al., 2020).However, in some situations, which tend to involve horizontal inhomogeneities, γ 0 > 45 • ; these include e.g., baroclinity, terrain-induced turning (especially with stability), convective cells, and various persistent storm structures. 6In the Southern hemisphere, the signs are reversed: geostrophic flow around a local low-pressure moves clockwise, with surface-induced turbulence ('friction') causing the flow to again increasingly turn towards low pressure as the surface is approached, and thus γ < 0.
The more generic form of this, for an arbitrary coordinate system, follows from (11): We note that ( 14) and ( 15) are more direct alternatives to dealing with functions of φ G , which become apparent if one expands cos γ in ( 12) or ( 13).However, one can see that there can be an angular dependence within the stress-related parts written above; when considered in coordinates defined with the x-direction aligned with the mean wind at height z, in the general forms ( 12) and ( 15), U/|S| and V /|S| can be written as cos φ and sin φ, respectively.From (12) and using cos γ = cos φ cos φ G + sin φ sin φ G , then in coordinates again defined by |S(z)| = U (z), after some rearranging we arrive at an expression for veer like ( 14): Compared to (14) this lacks a negative sign, but sin φ G is negative and with larger magnitude than the positive contribution to the denominator, ∂⟨uw⟩ ∥ /∂z/(f |G|); this will become more apparent in the sections which follow.We also remind that in these coordinates φ G = γ G (z), and opposite signs will occur for the southern hemisphere (expressions 14-16 give dφ/dz signed for the northern hemisphere in mathematical coordinates, reflecting winds rotating on average clockwise with increasing height, i.e. negative veer).
For wind energy ∂(cos φ)/∂z might be considered as relevant as ∂φ/∂z, because it allows direct expression of the veerinduced variation in streamwise wind velocity component relative to a reference height such as hub height.One could expect that the reduction of cos φ away from a given z counteracts the effect of typically positive shear; if desired, the veer can be simply re-expressed later in terms of cos φ for a given coordinate system, instead of trying to use an expression such as (12).

Misalignment of shear and stress
One can see a connection between the shear, veer, and stress in (9) and ( 12), and we can further examine the relation between shear and stress using complex notation as in (7).The 'misalignment' can be expressed via the angle between ∂S/∂z and ⟨sw⟩, i.e.The root of such misalignment arises in the rate-equation for ⟨sw⟩.In the limit of horizontal homogeneity, if we combine the stress budgets (e.g.see Wyngaard, 2010), i.e. adding ∂⟨uw⟩/∂t to i∂⟨vw⟩/∂t, we may write The pressure-strain contribution has been written as ⟨sw⟩/τ R via the commonly-used Rotta parameterization, where τ R is the Rotta time scale; this is the basis for commonly used flux-gradient relations (Wyngaard, 2004).In such mixing-length relations, i.e. using the 'Boussinesq hypothesis,' ⟨w 2 ⟩τ R is simply written as a turbulent diffusivity −ν T , and the final term in ( 18) is neglected.We continue to neglect advection and horizontal transport (such as U ∂⟨sw⟩/∂x and ∂⟨suw⟩/∂x respectively); these can also contribute to misalignment between ∂S/∂z and ⟨sw⟩ in areas of upwind horizontal inhomogeneity such as nonuniform terrain and turbine wakes.Thus in models where an eddy-diffusivity (flux-gradient relation) is used, such as most RANS solvers which employ 2-equation turbulence models, for flow over homogeneous surfaces there will be no stress-shear misalignment.
Ghannam and Bou-Zeid (2021) derived a dimensionless relation in terms of the angular differences β ma and γ instead of velocity components; although it does not afford convenient description of the veer, it can be re-cast to show the effect of the misalignment angle: Thus when the stress is aligned with the shear (β ma = 0), then f |G| sin γ = −∂|⟨sw⟩|/∂z; this can be seen as a case of (13).
The contribution of stress-shear misalignment to the veer can also be seen considering ( 18) with our earlier derivations, with misalignment modifying the stresses.For example the cross-wind stress in ( 13)-( 15) can be written since the Rotta timescale can be expressed in terms of turbulent kinetic energy k via ν T = τ R ⟨uu + vv + ww⟩/3 (see Pope, 2000;Hatlee and Wyngaard, 2007).But the turbulent third-order moment ⟨sww⟩ is difficult to measure, so a model for it would be needed in order to explicitly incorporate misalignment into veer predictions.Fortunately the misalignment β ma tends to be small in the surface layer (Geernaert, 1988), and also beyond the surface layer over homogeneous terrain or long fetch over water, especially without baroclinity (Berg et al., 2013).However, it has been known for decades (Moeng and Wyngaard, 1989) that turbulent transport is relevant in convective ABLs, so one expects more misalignment in unstable conditions; indeed Santos et al. (2021) saw this from measurements over multiple heights over a land and sea site, as did Berg et al. (2013) to a lesser extent (due to the relatively short measurement campaign) over water.The misalignment tends to be smaller in neutral conditions, and thus we do not (yet) offer explicit treatment for it.

Alignment of shear and stress; canonical solutions
When turbulent transport of stress is negligible (along with baroclinity and inhomogeneity), in steady conditions the stress and mean velocity gradient are aligned, allowing use of an eddy diffusivity.The veer can then be cast as a nonlinear differential Ekman solution.
Ekman (1905) assumed the turbulent stress was related to the mean shear using a constant eddy-viscosity ν Ek , which in our notation is expressible as ⟨sw⟩ = −ν Ek ∂S/∂z.Thus the momentum balance (7) simplifies to which gives the classic Ekman solution where the characteristic Ekman (e-folding) height h Ek is defined as h Ek ≡ 2ν Ek /f .Simpler than relating Ekman veer to shear, the solutions above along with (9) give the veer directly as  (24) this result has units of radians/m measured counter-clockwise, with the linear approximation7 deviating from the exact form by less than 1% for z < 1.5h Ek .Integrated over z± ∆z/2, this gives the veer across an extent ∆z: The Ekman forms might be seen as an upper limit on veer for h Ek on the order of typical ABL depths (∼300-1000 m), analogous to what was found by van der Laan et al. ( 2020) for the cross-isobar angle γ 0 .
From (23) one can also find an expression for the Ekman shear exponent α Ek via (1), This may also be seen as an upper limit, particularly in the surface layer where an unrealistically large diffusivity is assumed; one can see that Ekman theory predicts α → 1 approaching the surface.
Linear diffusivity profile: 'modified Ekman' or surface-layer regime Using a surface-layer eddy-viscosity relation ν T (z) = κu * z consistent with ASL theory, Ellison (1956) derived the solution of ( 7), resulting in a profile of wind vector 'deficit' expressible as (Krishna, 1980) where ker 0 (x) and kei 0 (x) are the so-called Kelvin functions (see e.g.Abramowitz and Stegun, 1972).But the Ellison solution can be written more compactly and conveniently, similar to ( 23) with a complex argument, as K 0 (x) is the zeroth-order modified Bessel function of the second kind, and the modified-Ekman length scale is defined by For the range 0.02 ≲ c G ≲ 0.06 encountered in nature under neutral conditions (Hess and Garratt, 2002;van der Laan et al., 2020), for z H /h mE ≫ 0.1 the arctangent of ℑ{S}/ℜ{S} can be approximated via series expansions of ( 27) or ( 28) to yield the practical result this follows the numerical solution to within ∼20% for 0.3 ≲ z/h mE ≲ 2, moreso for c G approaching 0.04.
It was shown in van der Laan et al. ( 2020) that the Ekman and Ellison solutions basically gave upper and lower limits, respectively, to observed full-ABL turning (φ G − φ 0 ).Following this, in Fig. 1 we present veer profiles along with the relationship between veer and shear, for the Ekman and Ellison solutions; the former is calculated via the expressions in ( 24) and ( 26), while the latter is obtained via (28).One can see in Fig. 1 that the Ekman solution produces effectively less mixing away from the surface (for z > ν Ek /(κu * )), and consequently a higher shear exponent than Ellison's.Similarly, the dimensionless Ekman veer exceeds that predicted by the Ellison solution, for z/h Ek ≳ 0.1, consistent with γ 0,Ek = 45 • , which exceeds the γ 0 of 5-15 • predicted by Ellison's (van der Laan et al., 2020).However, we note that the depth h can differ for Ekman and Ellison solutions; h Ek = h mE only if one chooses ν Ek = (κu * ) 2 /2f .We also point out that for larger c G , i.e. smaller Ro 0 (not shown in figure), near the surface (z/h mE ≲ 0.1) Ellison's veer grows yet larger than the peak value shown at z ≈ 1.5h mE and relative to the behavior seen for c G = 0.04; however, this idealized near-surface behavior is likely not relevant for wind applications.

Practical forms and application
To use the expressions derived for veer earlier, one needs the vertical derivatives of stress (or its profile) and the geostrophic wind speed; in particular the first and second vertical derivatives of the cross-wind stress ⟨vw⟩ ⊥ appear in ( 14) and ( 16), along with |G|.In wind energy applications, engineers typically lack site-specific stress profiles, unless they are taken from flow modelling; if the latter is reliable, then there is probably less need for the shear-based estimates for veer given in this work.
The large-scale horizontal pressure gradients which drive ABL flow, expressible as the geostrophic wind G, are likewise rarely measured (though lidar measurements above the ABL can make this possible, e.g.Pedersen et al., 2013).The shear contribution to veer is multiplied by |S/G| in ( 14)-( 16).To obtain a practical form relating shear to veer, we can start by parameterizing |S|/|G|; fortunately |G| is commonly calculated in practice using a geostrophic drag law ('GDL', Rossby and Montgomery, 1935).Long used in wind applications such as WAsP (Troen and Petersen, 1989) and related wind resource software, it is expressible in scalar form as where the empirical coefficients {A, B} are assumed to be constants in typical wind application.The geostrophic drag coefficient c G ≡ u * /|G| and ABL turning (cross-isobar angle) φ G are seen to vary with surface-Rossby number Ro 0 (Blackadar and Tennekes, 1968); these and {A, B} have been shown to depend on dimensionless stability L −1 u * /f (Arya, 1978;Kelly and Troen, 2016), strength of ABL-capping inversion (Zilitinkevich and Esau, 2002), and baroclinity (Arya and Wyngaard, 1975;Nieuwstadt, 1984).For practicality, we start by assuming near-neutral stability, which is appropriate in the mean for most places, as it represents by far the most frequently observed conditions (Kelly and Gryning, 2010); we continue to neglect baroclinity; and we neglect influence of the capping inversion strength.(approximate) 'reverse' form of (30) to get the drag coefficient as (Troen and Petersen, 1989) where the surface Rossby number is Ro 0 ≡ |G|/(f z 0 ) and c rGDL is taken to be 0.485 following its use in the wind resource program WAsP for several decades.Alternate forms of ( 32) exist, such as that of Hess and Garratt (2002); the latter corresponds simply to setting A = 1.28 and c rGDL = 0.472 in (32).For a given roughness length z 0 and measured wind speed |S|, lacking the (surface) friction velocity u * , one needs a relation to connect u * and |S|, in order to get |G|.This can be done through the same wind profile relation upon which the GDL is built, i.e. the log-law; one can use u * = κ|S|/ ln(z/z 0 ) within ( 30) or alternately |S|/|G| = (c G /κ) ln(z/z 0 ) using ( 32), where in the latter ( 30) is also employed to find |G| within Ro 0 .
In practice one would like a direct estimate for the veer, using the routinely-measured shear, since α is seen to drive ∂φ/∂z.
One way could be to just ignore the stress divergence terms in ( 14) or ( 16), which with calculation of |G| mentioned just above considerably simplifies the problem.However, this might not be justified, particularly if u 2 * /(f h) is not negligible compared to |S|, as seen from comparing contributions to ( 14)-( 16); this can be seen using the scaling ∂⟨uw⟩/∂z ≈ u 2 * /h where h is the ABL depth (e.g.Wyngaard, 2010).Thus we consider estimating vertical derivatives of the stresses, starting with the ∂⟨uw⟩/∂z just mentioned, which can be used in ( 16).Similarly, one can estimate where Ro h ≡ G/(f h) is the Rossby number based on ABL depth and c vw is of order 1; we will treat c vw as an empirical constant which is tuned later below.To use ( 16) we also need to find sin φ G ; employing (31) and using trigonometric identities to expand sin(φ − φ 0 ), with some rearrangement one obtains Employing this, (33), and ∂⟨uw⟩/∂z ≈ u 2 * /h, along with with (30) or (32), allows one to then use ( 16).On the other hand, using ( 14) is simpler and more convenient than ( 16), because it only requires ∂⟨vw⟩/∂z in addition to the second derivative of ⟨vw⟩ just approximated in (33) above, so one can also simply approximate ∂⟨vw⟩/∂z ≈ h∂ 2 ⟨vw⟩/∂z 2 and use (33); the GDL forms ( 30) and (32) then allow one to get Ro h and c G , respectively.Whether using ( 16) or ( 14), we note that the shear contribution to veer includes a surface-Rossby number (Ro 0 ) dependence through S/|G|, while the stress-Coriolis contribution includes an ABL-depth dependence, Ro h ; either way, if we do not neglect the latter, then we also need an estimate for the ABL depth h.If the shear contribution is expected to dominate variations in veer, then the estimate of h may not be so crucial; we will consider this further below in our comparison with real-world cases, and also direct interested readers to e.g.Liu and Liang (2010) for statistics of h in different conditions.This section presents analysis of results from RANS simulations of the neutral atmospheric boundary layer9 , and of observations at different sites (which include the impacts of stability).The simulations are analyzed to check the relations given here, as well as examine the behavior of and contributions to veer across the range of Rossby numbers (Ro 0 and Ro h ) encountered in nature.Investigation of observations, spanning turbine rotor heights in five places having different wind regimes and conditions, includes probing the interconnected behaviors of shear (exponent) and veer with atmospheric stability -as well as their joint statistics, universal trends, and variation with wind speed.The statistical demonstration of observations is accompanied by predictions of veer using empirically updated forms of the relations given in the previous section, as well as the forms themselves.

Model and setup
The Navier-Stokes solver Ellipsys1D (van der Laan and Sørensen, 2017), which is a one-dimensional version of the multiblock general CFD solver Ellipsys3D (Sørensen, 1995), was used to simulate the Reynolds-averaged flow in neutral atmospheric boundary layers, including Coriolis forces.Assuming zero vertical velocity and constant pressure gradients, it solves the RANS equations for incompressible flow with a finite-volume scheme.The ABL 'top' (above which turbulence is extinguished) is modelled via the length-scale limiter model of Apsley and Castro (1997)  The domain height is set to 10 5 m to ensure it is much larger than h for all simulations, and the bottom boundary is handled by a rough-wall condition (Sørensen et al., 2007).The numerical 'grid' is a vertical line, with the bottom cell height being 1 cm (placed above the roughness length) and the cells' sizes growing progressively upward with an expansion ratio of 1.2; the total number of cells is 384.At the bottom cell a Neumann condition is set for k (dk/dz = 0) and ε is set to the logarithmic value, the wall stress is consequently defined by the neutral surface layer for this cell.More details, including a grid-refinement study, may be found in van der Laan et al. (2020).
Using a constant geostrophic wind speed, the flow is driven by a constant pressure gradient, starting with an initial wind profile set to |G| at all heights; the ABL depth grows upward until convergence occurs, providing a steady solution and h for a given choice of z 0 , pressure gradient (thus G and f ), and turbulence (k-ε) limiting lengthscale ℓ max .The Buckingham Pi theorem can be used to reduce the four parameters {z 0 , G, f, ℓ max } into two dimensionless groups, namely Rossby numbers for

Shear and veer over neutral ABLs simulated over entire range of Rossby numbers found in nature
First we check that the RANS simulations confirm the shear-veer relations developed earlier; we expect this to be, since there are no extra terms in the simulated Navier-Stokes equations compared to (7). Figure 2 displays both sides of ( 9) and ( 13) respectively, for four cases representing somewhat common real-world conditions, for heights between 50-200 m.From Fig. 2 • Figure 2. Demonstration of 1D RANS solver results, conforming to eq. 13 (left) and eq. 9 (right).Dashed line represents 1:1 prediction; simulated ABL depth heff calculated from (35).
one can see that the Ellipsys1D solutions conform to equations ( 9) and (13) derived earlier.
Towards investigating the behavior of veer (and shear) in terms of Rossby numbers -which is facilitated by RANS, but is quite difficult to accomplish with measurements -we turn our attention to the variation in veer as a function of surface roughness.Admiting that we are using one-dimensional simulations over a homogeneous surface, we now consider the directional change across typical turbine rotor heights, i.e. ∆φ from z = 50 m to 150 m. the right-hand plot in Fig. 3 one can see that ∆φ is roughly proportional to 1/ ln(Ro 0 ), as expected from the S/|G| contribution to veer considering ( 30) and ( 32).
Looking back on (33), we may also expect a Ro h dependence in the veer, due to the stress gradient contributions.only for one roughness, because the curves of veer versus ABL depth and Ro h look nearly identical when using any other z 0 (or Ro 0 ) value, such as water roughnesses less than 0.3 mm.In other words, the sensitivity of veer to h is essentially independent of z 0 , if one varies these separately from case to case as in our numerical simulations.Looking at these results, we note a behavior that is not inconsistent with the estimates for stress-gradient contributions following (33): the veer is proportional to Ro 1.4  h (or h −1.4 ) over a range of ABL depths routinely observed in reality (h ∼200-800 m), though the dependence softens to which are more rarely encountered, the height dependence vanishes; this can be intuitively interpreted, as ∆z/h becomes so small that less directional change is found for a given ∆z when h is increased further.The veer and its h-dependence is seen to be basically independent of height for these 50-100 m vertical spans: at the heights of interest for wind energy shown, the lines collapse onto one another.To compare with Fig. 3, multiplying the veers in Fig. 4 by ∆z = 100 m for the blue and gold curves, we can also see that for a realistic range of ABL depths and roughnesses, the effect of h is stronger than that of z 0 : across all Ro 0 a variation in ∆φ| 150 m 50 m of only several degrees is seen, whereas across the common range of Ro h a variation of more than 15 • is shown.Now that we have seen in Figs.3-4 how the veer (or simply the turning ∆φ for typical rotor ∆z) depends on z 0 and h, presumably due to the S/|G| (shear) and stress-gradient contributions respectively, it is prudent to examine the relative sizes of each of these contributions -particularly because RANS affords us this opportunity.One can cleanly separate these contributions by examining the variation of cos γ, as indicated by ( 12) and ( 13).Accordingly, Fig. 5 presents the two contributions to the dimensionless veer ∂ cos γ/∂z derived in ( 12), for the four over-land cases shown in Fig. 2 as well as an over sea case with the same ABL depth as two of the land cases.One can note from Fig. 5 that the shear and stress-gradient/Coriolis contributions largely offset each other, with each being an order of magnitude larger than their sum, which is equal to the dimensionless veer ∂ cos γ/∂z.The vertical profiles of 'pointwise' veer shown in the figure, which were calculated using 3rd-order finite difference, indicate that in neutral conditions the veer is smaller offshore compared to on land.Further, one sees the combined effect of the behaviors noted from the previous two figures: shallower ABLs have larger veer, as do ABLs over rougher surfaces, with Ro 0 (z 0 ) having a smaller impact than Ro h (h).
This can be put into a more practical context by considering the variation of shear and veer together across the range of Rossby numbers found in atmospheric flows.From the left and center panels of Fig. 6, it becomes evident that Ro h affects ∆φ more than α for typical rotor extents; opposite of this, Ro 0 affects the shear more than the veer.Further, for the relatively representative set of (common) cases shown in the center and left-hand plots in Fig. 6, we notice much less variation in {α, ∆φ} compared to the entire parameter space displayed in the right-hand plot; as we will see in the next sub-section, the right-hand plot is more in line with observations, despite the RANS solutions representing nominally neutral conditions 10 over uniform surfaces with neglect of shear-stress misalignment and baroclinity.

Results from measurements in different wind regimes and sites
After examining the behavior of neutral-ABL dependencies for shear and veer above from simulations, now we consider the behavior of each in the real world from measurements at different sites, which includes e.g. the affects of stability.The datasets are the same analyzed by Kelly et al. (2014a), which showed shear exponent statistics for these locations, except a longer record of Høvsøre data was used for the current study (10 years, from 2005-2015).These are: the aforementioned Høvsøre site, from 60-160 m height for both homogeneous land and sea sectors; the partly forested but flat Østerild site (Hansen et al., 2014) for two virtual rotor spans, from 45-140 m and 80-200 m over one year; the Dutch research site Cabauw (Beljaars and Bosveld,   10 One could argue that our RANS solutions can also be interpreted to include stable conditions, since the lengthscale-limited k-ε turbulence model can have its maximum mixing length ℓmax rewritten using the Blackadar (1962) mixing-length formulation such that ℓ −1 max,eff = ℓ −1 max plus a stability contribution, as shown in van der Laan et al. (2020).However, such interpretation employs M-O theory along with the Blackadar-type form to 'combine' a surface-layer scale ℓ ASL ∝ z with ℓmax; here we choose to keep our analysis as general as possible -avoiding particular ASL forms or assumptions, and models for turbulence length scale.We investigate the statistical behavior of veer with shear exponent as well, not only to see their interdependent behavior, but also towards providing useful relations for their variability and practical prediction of veer from typical wind energy measurement campaigns.

Shear exponent
Here we briefly explore the connection between probability distribution functions (PDFs) of stability and shear exponent.The shear distribution f (α) can be connected to f (L −1 ) in the surface layer during stable conditions, but there is not necessarily a one-to-one (unique) mapping between the two (Kelly et al., 2014a).As seen in ( 4) and ( 5), α tends to correlate with stability (1/L) and particularly buoyant destruction (−B) during stable conditions, when turbulent transport is negligible.This is shown in Fig. 7, which displays the joint probability density of α| 160 m 60 m and L −1 calculated in the ASL at z = 10 m from the homogeneous land sectors at the Danish national test station of Høvsøre (Peña et al., 2016) from 10-minute averages over a 10-year period.), stable (L −1 > 0.001m −1 ), and unstable (L −1 < −0.001m −1 ) flow regimes, weighted by frequency of occurrence to show the relative contributions to the overall distribution.The threshold of ±0.001m −1 for L −1 is a sensible choice because then z/|L| ≪ 1 (consistent with neutral conditions) in the surface layer, which is generally taken to have a thickness of 100 m or less (roughly h/10, also recalling that M-O theory's applicability diminishes with height above the surface-layer).Even at this relatively flat and uniform site, negative shear happens in both stable and unstable conditions, though moreso in unstable and yet less often in neutral conditions; overall, α < 0 occurs less than 5% of the time over 60-160 m here, and 8-9% from anemometers at 100-160 m heights (not shown).We also note that while the 'ideal' Høvsøre land (eastern) sectors have conditions split somewhat evenly between the three stability regimes, other sites can differ (Kelly and Gryning, 2010).

Veer
Along with distributions of α, measured veer distributions are shown in Fig. 8 for both land and sea conditions at Høvsøre, i.e.From the two plots one can see that the most common α and ∆φ/∆z, i.e. the portions of P (α) and P (∆φ/∆z) with respective probabilities within an order of magnitude of the peak values, both systematically differ when using higher measurements at 100-160 m compared with 60-160 m heights; however, the shift in the commonest α is significantly smaller than the analogous shift in ∆φ/∆z between these two height ranges.This happens over both land and sea, though both α and ∆φ/∆z 100-160 m is less than +5% over land and < +30% is seen over sea, while the mean veer ⟨∆φ/∆z⟩ is seen to increase by factors of ∼ 5/3 and 2 over land and sea, respectively.
There are several other notable differences between the shear and veer statistics shown in Fig. 8.The peak portion of P (α) is significantly wider over land compared to offshore (with larger σ α over the rougher surface), while the shape around the P (∆φ/∆z) peak does not differ significantly from land to sea here.Further, the (logarithmic) slope of P (∆φ/∆z) versus ∆φ/∆z for veer larger than the PDF peak is basically the same regardless of height or surface conditions; this and the land-sea difference between P (α) are consistent with the earlier RANS results, where z 0 primarily affects α, while ∆φ/∆z is impacted more by ABL depth.P (α) also has wider 'tails' (extremes) higher from the ground on both sides, including negative shear due to low-level jets (such as that due to the capping-inversion when h ∼ 200 m), whereas the veer simply becomes larger due to such jets in shallow ABLs, as jets and the environment associated with the capping inversion simply causes more turning, and not a reversal.The negative veer occurs due to nonstationary processes like passing fronts (e.g.Clark, 2013), as well as baroclinity and motions associated with it (Arya, 1978;Foster and Levy, 1998;Floors et al., 2015).Comparing the solid and dashed lines in Fig. 8, one sees that the highest veers ∆φ/∆z are larger for the 100-160 m measurements than those from 60-160 m; this is again due to more impact of the ABL-capping inversion and associated jet with turning.
As with shear, stability affects veer, with stable conditions expected to lead to higher veer due to its damping effect on vertical fluxes (suppressing vertical 'communication' of flow information).Following the plots shown in Fig. 7 for the shear exponent α, Fig. 9 displays the effect of stability on veer for the Høvsøre land sectors.The figure shows P (∆φ/∆z) for neutral (|L −1 | < 0.001m −1 ), stable (L −1 > 0.001m −1 ), and unstable (L −1 < −0.001m −1 ) flow regimes, weighted by frequency of occurrence (indicating relative contributions to the full PDF), as well as the joint distribution of stability and ∆φ/∆z.From Fig. 9 one sees that in comparison with P (α) shown in Fig. 7, the peaks of veer distributions P (∆φ/∆z) do not depend so much on stability.However, as with the shear distribution, ∆φ/∆z also has its largest values dominated by stable conditions; this makes sense considering that stability tends to maintain vertical gradients by limiting vertical fluxes.Unlike the results shown for the RANS simulations or predicted by theory, negative veer occurs as in Fig. 8 and described thereunder; one can see in Fig. 9 that it basically happens during non-neutral conditions, which tend to occur at lower wind speeds, and is dominated by unstable conditions.Looking at the joint distribution P (α, ∆φ/∆z) one sees that for the most common veer values (0 ≲ ∆φ/∆z ≲ 0.1 • /m), which tend to occur around neutral conditions, there is a mild stability dependence; however for less neutral conditions there is little correlation between veer and stability, aside from higher veer simply being observed more often in stable conditions.

All
To show the behavior of veer across different locations, Fig. 10  most commonly observed values (which tend to occur in stable conditions, as shown in Fig. 9 for the Høvsøre case above), the distributions behave similarly across locations; in particular the "slope" of the semi-log plot for veer exceeding the PDF peaks is roughly constant for ∆φ/∆z ≳ 0.2 • m −1 in each case.These slopes correspond to (conditional) PDFs for the largest veer of the form where the charcteristic veer scale defined by Υ −1 veer ≡ ∂[ln P (∆φ/∆z)]/∂(∆φ/∆z) ranges from roughly 0.07 to 0.11 • •m −1 .The lowest Υ veer corresponds the offshore Høvsøre case, while the highest Υ veer matches the Østerild case from 45-140m.We expect larger Υ veer to correspond to occurences of higher 1/L, i.e. a larger width σ + of the stable-side distribution P (1/L) following Kelly and Gryning (2010); essentially the large-veer PDF in ( 36) is conditional on stable conditions, i.e. we could express it as P ∆φ/∆z L −1 > 0 ∝ exp [−(∆φ/∆z)/Υ veer ].The dominance of stable conditions reported by Peña (2019) https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.
for z ≳100m at Østerild is consistent with this, though the data from z =80-200m (green line in Fig. 10) with smaller apparent Υ veer might appear to not be, considering the increasingly stable conditions higher up at this site; but looking at the Østerild curves in the figure we see that for higher veer ∆φ/∆z ≳ 0.4 • •m −1 , there is consistency: the two largest Υ veer occur for z =80-200m and z =45-140m, respectively.Future work needs to be done to explore this, since we lack air-sea temperature differences (or water-air heat flux) for the Høvsøre offshore case and stability information for the MR site, while stability effects above forests tend to be diminished and are difficult to interpret due to turbulent transport through the treetops (e.g.Sogachev and Kelly, 2016).
One also sees the peaks of P (∆φ/∆z) in Fig. 10 are at smaller ∆φ/∆z for the forest-dominated Østerild cases, with the peak of the offshore Høvsøre veer distribution falling between these and the ∆φ/∆z corresponding to the the land cases of Høvsøre, Cabauw, and 'MR'.We remind that the most commonly-found veer values are generally dominated by neutral conditions (or modestly stable for the exceptional Østerild site above 100m), and point out that the mode of ∆φ/∆z is essentially the same (0.005-0.006 • • m −1 ) for the land cases that are not dominated by forest.Further considering the RANS simulation results from Fig. 6 discussed earlier, the mode of ∆φ/∆z being smaller for Høvsøre offshore than for the land cases (of Høvsøre, Cabauw, and 'MR') can be explained by the smaller ABL depths most commonly observed offshore compared to onshore; this is consistent with the ABL depth distributions aggregated and reported by Liu and Liang (2010).The mode of ∆φ/∆z found at the inhomogeneous forest-dominated site Østerild are more strongly affected by the tree-enhanced mixing (which reduces the veer magnitudes) and to a lesser extent by shallower ABLs due to the coastline 5-20km upwind in some directions.
The dependence of veer on wind speed at the sites considered is shown in Fig. 11, which displays the joint distribution of veer and 10-minute mean wind speeds, P (∆φ/∆z, U ). Along with the joint distribution, the mean veer conditioned on wind speed, ⟨∆φ/∆z⟩ U , is displayed.
From Fig. 11 one can see results consistent with the effects of stability discussed earlier and evoked by Fig. 9: at higher speeds neutral conditions dominate, giving decreased mean veer.This is more pronounced for the onshore cases (though there is still a reduction of nearly 40% going from 12 to 24 m s −1 for the offshore case), because sea-air heat fluxes and associated 1/L magnitudes tend to be relatively smaller due to water's large heat capacity (e.g.Cronin et al., 2019).It is notable that for the representative wind turbine rotor heights considered, the veer tends to be largest for wind speeds below typical turbine rated speeds, especially over land; this can have consequences on both the power output and effective power curve for preconstruction AEP estimates, as well as loads.
Further, a narrower range of veer with increasing wind speed is seen in Fig. 11, regardless of surface properties; such narrowing is impacted by stability, but also occurs in neutral conditions.The variability of veer with mean wind speed is presented in Fig. 12, which displays the standard deviation of veer conditioned on mean wind speed for the sites/cases considered.It also adds a line to show the overland Høvsøre case filtered for neutral conditions.
Consistent with the joint-PDFs P (∆φ/∆z, U ) in Fig. 11, from the semi-logarithmic plot of standard deviation of veer conditioned on mean wind speed in Fig. 12 we can see that the variation in veer decreases with wind speed, and moreso over land than water.It is also seen that for the onshore Høvsøre case σ (∆φ/∆z)|U is smaller in neutral conditions compared to over all stabilities, with the two values converging at higher speeds due to the increasingly neutral conditions.For each site having a standard deviation of veer over all speeds σ ∆φ/∆z and mean wind speed ⟨U ⟩, the rms veer conditioned on wind speed roughly follows the empirical form up to about 12 m s −1 over land, and to higher speeds offshore.A more complicated speed-dependent variability in veer is 530 seen for the MR case, with higher σ (∆φ/∆z)|U at speeds above 15 m s −1 caused (presumably) by hill-induced turning.This has two consequences worth mentioning: first, that turbines at a site such as 'MR' can experience persistent veer above rated speed, potentially increasing loads and/or reducing power below rated; secondly, such speed-dependent behavior is likely difficult to capture with standard single RANS simulations, demanding more detailed treatment to handle the Reynolds-number dependence despite the lack of stability effects at such speeds.

Relating veer to shear in application
One of the aims of this work is to relate veer to shear (or shear exponent), as with the expressions developed in Sect.2.3-2.4.
Here we present joint observations of shear exponent and veer, and following these, give practical simplified forms based on the equations derived earlier in sections 2.3-2.4.
Following the previous subsection, we first consider the joint behavior of ∆φ/∆z and α with wind speed and stability, for the 'simple' onshore Høvsøre case having homogeneous upwind conditions.Figure 13 shows the observed joint distribution P (∆φ/∆z, α) in neutral conditions, over typical turbine operation speeds (4-25 m s −1 ) and separately over different speed ranges (4-8, 8-12, 12-16, and 16-25 m s −1 ); counts are used instead of PDF per wind speed range, to show relative frequencies of occurence.
From Fig. 13 one can notice that in neutral conditions there does not appear to be significant variation in the joint shear-veer behavior with U , with a bit more variability at the lowest speeds and smaller values of both ∆φ/∆z and α for U > 16 m s −1 ; this is consistent with Figs. 7, 9, 11, and 12.The larger spread at lower speeds for neutral conditions is attributed to the larger relative effect of non-stationarity and particularly sampling uncertainty; per the latter the integral time scale increases roughly as U −1 (Wyngaard, 2010) so fewer integral time scales are 'sampled' per each 10-minute period.This is also evident considering the previous plot of σ (∆φ/∆z)|U versus U in Fig. 12, where one sees σ (∆φ/∆z)|U increasing with diminishing wind speed during both neutral and all conditions for the Høvsøre land case, but where stability effects cause larger veer variability up to speeds of about 15 m s −1 .Also, the overall jPDF P (∆φ/∆z, α) appears similar to that in the most common speed range (8-12 m s −1 ).Aside from nonstationarity and sampling effects one does not expect much speed dependence in neutral conditions, considering the α-related part of ( 14)-( 16) behaves as |S|/G, which following (32) has a weak |S|-dependence through (ln Ro 0 − A) −1 ; the RANS results also confirm this.We note a joint trend between α and ∆φ/∆z, but also see a spread around the most common shear exponent and veer values due to variations in ABL depth, stress gradient and curvature, and top-down stability (capping-inversion strength, see e.g. Kelly et al., 2019a), in addition to nonstationarity.Figure 14 shows joint α-veer distributions like Fig. 13, but over all conditions, i.e., not limited to neutral stratification.One notices immediately the more frequent occurence of higher veer and shear, as well as negative α and ∆φ/∆z.Further, in addition to a wider range of shear and veer compared to neutral conditions, in Fig. 14 one can see there is also a sharper increase in ∆φ/∆z with α for larger α, due to stable conditions.One can see that at the most common (8-12 m s −1 ) and lower wind speeds, which occur in the range below rated speed for typical turbines, there is a significant increase in ∆φ/∆z with α in the more stable conditions where α ≳ 0.3; this higher 'slope' of ∆φ/∆z vs. α is likely enhanced by the shallower ABLs which generally occur along with stable surface-layer conditions (we remind that the stability metric L −1 was measured in the ASL), whereby additionally stable air above augments the veer.As mentioned previously, the turning and veer near the ABL top will continue to increases for yet shallower ABLs (decreasing h); meanwhile α is less sensitive to h as the upper height (used to calculate α and ∆φ) exceeds the peak of the inversion-induced 'jet.' Further, such high-veer conditions are not rare for such a 'simple' site at the heights considered (60-160 m); e.g., conditions where ∆φ/∆z = 0.2 and α = 0.4 (a veer of 20 • over a 100 m rotor) occur as frequently as conditions with zero shear and veer.
Towards relating veer to shear for application, we now consider the mutual behavior of ∆φ and α together at all of the sites analyzed for this work.Figure 15 shows the joint distribution of shear and veer for the sites considered, with each plot also including the conditional mean of veer per shear exponent (i.e.⟨∆φ/∆z|α⟩, as solid lines).From this figure we see a number of trends across the six cases analyzed.First, some nonlinear variation of veer with α is evident, along with the (less common) occurence of negative values of shear and veer, as was seen in Fig. 14 for the Høvsøre land case.Further, the veer tends to be skewed towards higher values: i.e., ⟨∆φ/∆z|α⟩ exceeds the most commonly observed values of ∆φ/∆z; however, the site MR does not show such skewed behavior (consistent with Fig. 10), presumably due to the complex terrain there.We note the conditional mean veer ⟨∆φ/∆z|α⟩ is also more clearly nonlinear in α, becoming less dependent on α in low (and negative) shear conditions; the site MR is an exception to this, with hill-induced height-dependent turning causing larger veer for α smaller than the most commonly-observed values there.

Simplified estimate of veer per α
Figure 15 also includes two predictions based on the theory presented earlier.First, as discussed at the end of section 2.4, using only the shear-associated (|S|/|G|) portion of ( 14) to be practical, we arrive at the estimate Lightest shades are 10% as likely as darkest shade in each plot.
above the most common α observed for the onshore sites considered (α ≳ 0.2) and at α ≳ 0.1 for the offshore Høvsøre case; here we have used effective roughness lengths consistent with earlier studies employing these sites (z 0 =1.5 cm for Høvsøre land, 3 cm for Cabauw, 0.9 m for Østerild, 2 m for MR, and 0.02 cm for offshore).A value of c sα = 0.5 can be seen to fit the heterogeneous terrain cases where terrain and roughness dominated over stability (Østerild and MR, bottom plots of reduces the reverse drag-law constant by roughly 10-40% for the Rossby numbers applicable at these sites, consistent with the values of c sα used in the plots of Fig. 15; but again, to model stability effects beyond the surface layer becomes rather complicated and is the subject of ongoing work.For reference, a value of c sα = 0.6 fits the mean veer for the Høvsøre land case during neutral condtions (not shown), in contrast to the value of 0.7 which fits when all stabilities are considered there.

Veer estimate including both α and cross-wind stress
We remind that for simplicity, (38) ignored the effect of cross-wind stress; it neglects not only ⟨vw⟩ but consequently also Ro h , though it does incorporate the effect of Ro 0 seen in the simulations of Section 3.1.Thus we also consider an approximation of the ⟨vw⟩ terms using ( 33) in ( 14), which introduces Ro h , along with the parameterization for |S|/|G| from (38): where c G is found using (32), |S|/|G| is calculated the same way as done earlier for (38), and |G| within Ro h is calculated as it was within Ro 0 of (38).To use (39) the ABL depth must be prescribed, along with the constant c vw and the parameters {z, |S|, z 0 , f } also employed for (38).Given the negative curvature of lateral stress, ∂ 2 ⟨vw⟩/∂z 2 <0 (e.g.Wyngaard, 2010), c vw is negative and of order 1, with the ⟨vw⟩ (Ro h ) contribution reducing the predicted veer compared to (38).With its moderating effect on the α contribution, the ⟨vw⟩ part can produce an α-dependent 'upturn,' though slight; this is seen for the offshore and MR cases in Fig. 15.However, the constant c sα within |S|/|G| needs to be increased in order for (39) to fit the observed ⟨∆φ/∆z|α⟩; the values of 0.5 are replaced by 0.7, and c sα =0.7 and 0.8 for Høvsøre and Cabauw are replaced by 0.8 and 0.9, respectively.The value of c vw giving the estimates shown in Fig. 15 was −0.7 for all sites, while characteristic ABL depths h were taken to be 800 m over the simple land cases, 600 m offshore, and 1000 m over the hilly/forested terrain cases; we note that the results have limited sensitivity to h, but choose these values to be consistent with mean ABL depth observations over sites of similar character and h distributions aggregated by Liu and Liang (2010).One can see from Fig. 15 that the estimates of ⟨∂φ/∂z|α⟩ using (39) are not better than the simpler form (38), though the constants c sα and c vw could easily be 'tuned' together to give a better fit for each case.However, in practice one might not be able to do so, and wishes to simply predict veer based on α; to this end, for practical applicability we suggest using (38).Though such a recommendation would appear to be neglecting Ro h and the ABL depth, we note that for estimation of mean veer (per shear) one is not so concerned with variations of Ro h or Ro 0 at a given site.The spread (scatter) around the mean veer seen in Fig. 15 is due to variation of stability as well as Ro h or Ro 0 , and variation from site to site is also due to different distributions of Ro h or Ro 0 ; this is consistent with Fig. 6 and discussions following it.
To illustrate the differences just mentioned, both the mean and standard deviation (spread) of conditional veer is shown in Fig. 16, for all the sites/cases considered .One immediately sees the character of ⟨∆φ/∆z|α⟩ tends to follow the type of site; offshore has larger veer for high α, simpler sites like Cabauw and Høvsøre onshore exhibit modest veer for large α, and the more complex sites have more limited veer for α around or above its most common values of α.But we remind that Fig. 15 shows that high-shear conditions offshore are relatively rare, and that α exceeding ∼0.3 is more common at the complex sites.We also see from Fig. 16 that for low-shear conditions (α <∼ 0.1), the simpler sites exhibit higher mean veer than offshore and yet more compared to the forested cases, while much larger veer is present due to upwind hills at the MR site for such low shear conditions (though somewhat uncommon, as seen in Fig. 15).From the middle plot we further note a that the long-term variability in veer σ ∆φ/∆z|α is lower offshore for the most commonly occuring α there, while veer variability does not differ so much for the most common conditions across the other sites/cases-except for the 45-140 m (lower) height range at Østerild, which shows larger veer variability due to being in the roughness sublayer above the forest there.In very high-shear conditions (α >∼ 0.5) the veer variability is highest offshore (though rarer).However, as shown in the right-hand plot of Fig. 16, the relative veer variability σ ∆φ/∆z|α /⟨∆φ/∆z|α⟩ tends to more clearly show the different character of the sites: the spread of veer relative to its mean (conditioned on α) is much larger in low-shear conditions over forest, while this relative spread is similar across all non-simple (forested, complex) cases for the most commonly occuring shear; the more homogeneous sites/cases exhibit comparable σ ∆φ/∆z|α /⟨∆φ/∆z|α⟩ under most conditions.For low-shear conditions, over more complex terrain the relative veer variability decreases, departing from the inhomogeneous forested (Østerild) values due to the large hill-induced mean veer.
The use of ⟨S|α⟩ in the calculations was also investigated; the plots in Figs.15-16 actually incorporated mean speed conditioned on α, though use of each site's corresponding overall mean speed ⟨|S|⟩ gave nearly identical results as those shown in the plots (within 2%, not shown).

Summary and Conclusions
We have derived relationships between shear exponent (α) and veer (∆φ/∆z), in a manner which avoids atmospheric surfacelayer (ASL) assumptions about meteorological parameters; this has been done in order to be applicable at wind turbine rotor heights, regardless of whether they are within or above the ASL.Canonical behavior of veer and shear with regards to surface The derived equations and RANS results essentially show that veer most simply arises from two contributions: the shear, and the vertical variation of crosswind shear stress at a given height (mostly through ∂ 2 ⟨vw⟩/∂z 2 , but also via ∂⟨vw⟩/∂z).
The numerical RANS solutions show that the shear and crosswind-stress contributions mostly offset each other in neutral conditions, and that each is much larger (up to an order of magnitude) than the veer itself.It is further seen that α primarily depends upon surface roughness in neutral conditions, with a weaker dependence on ∆z/h; in contrast, ∆φ/∆z more strongly depends on the ABL depth h, increasing as Ro n h where n is between 1 and 1.4 for the h most commonly encountered in nature (though ∆φ/∆z does also vary with 1/ ln Ro 0 ).These behaviors are consistent with the shear-veer relations derived in Sec.2.3.
We note that in this work we have also derived the cause of misalignment between shear and stress, as well as its contribution to veer; we remind that RANS solutions using mixing-length type closures (as well as e.g.WRF PBL schemes which lack turbulent transport) give stress aligned with shear, while the analytic shear-veer relations derived here allow for misalignment through the cross-wind stress.
The actual 'real-world' behavior of shear exponent and veer has also been investigated from multi-year measurements at four sites convering six different flow conditions (one with separate land and offshore sectors, one with measurements both in and above the roughness sublayer over a forest), for height spans or effective rotor diameters ranging from 47-60 m centered around (hub) heights of 88-140 m.The observed {α, ∆φ/∆z} include effects not fully accounted for in the equations derived here, particularly horizontal turbulent transport due to terrain inhomogeneities (Kelly, 2020) and nonstationary/transient flow conditions; though buoyancy is not explicitly accounted for, it primarily affects α and the stress, which are already incorporated into the derived veer equations.
The effect of surface-based atmospheric stability on shear and veer was examined for the relatively ideal (homogeneous) onshore site Høvsøre, where it is seen that unstable conditions dominate the low (negative) tails of the distributions P (α) and P (∆φ/∆z), while stable conditions are responsible for large α and ∆φ/∆z; neutral conditions contributed mostly to the peaks of the shear and veer distributions.Stability efffects are consequently seen to increase the long-term variability in veer and shear, as well as veer for a given α -particularly for the commonly-occuring wind speeds which tend to occur below the rated speed of modern wind turbines (e.g. Kelly and Jørgensen, 2017, Appendix B).The mean of both α and ∆φ/∆z was larger compared to neutral conditions, due to stably stratified conditions enhancing α and ∆φ/∆z more than unstable conditions (we note that sites having a distribution of 1/L more dominated by unstable conditions, possibly some offshore, could have mean behavior similar to that found in neutral conditions).
Comparison between offshore and homogeneous onshore sectors at Høvsøre showed α to be smaller offshore (as one would expect), with more extreme values at higher z (160 m) above the surface layer regardless of the surface; the latter is presumably due to the effect of the capping inversion for ABL depths which occasionally approach such heights (Liu and Liang, 2010;Kelly et al., 2014b).The veer distributions also show larger values over land compared to offshore, though to a lesser extent than P (α); but in contrast to α, which can increase or decrease (with wider extremes) due to the position of the jet associated Two practical veer-shear relationships were derived, including parameterizations for typically unmeasured quantities contained within them, then compared to the joint distributions P (α, ∆φ/∆z) and the ⟨∆φ/∆z|α⟩ measured from all sites over all conditions.A simplified form (38) neglecting the stress contributions was tested, as well as one (39) containing the cross-wind stress.Due to the relative simplicity of the practical shear-veer forms (and additional phenomena not included in them), they needed to be calibrated in order to match observed ⟨∆φ/∆z|α⟩; basically one coefficient in (38) and two in (39), all of which were of order 1 and universal (constant) across all six sites/flow situations analyzed.The form (39) for veer including crosswind stress did not give a better match to observations of ⟨∆φ/∆z|α⟩ across sites, compared to the simpler formula (38), and so we recommend the latter for shear-based predictions of veer at this time.Both forms provide their best predictions (within 10% of observed) ⟨∆φ/∆z|α⟩ during the most commonly-observed (moderate speeds and shear) and highest-impact (largeveer stable) situations, with underpredictions of mean veer occuring in low-veer conditions.The observed ⟨∆φ/∆z|α⟩ are nonlinear in α, whereas the derived forms were nearly linear, with the inclusion of cross-wind stress containing only a slight implicit nonlinearity.Lacking turbulent transport, our predictive mean veer relations are more suited for neutral and stable conditions where transport is less significant (e.g.Wyngaard, 2010); the underpreictions for smaller α, dominated by unstable conditions, evoke such.Consistent with this, the hilly site 'MR' shows yet more low-shear deviation from our predictions due to inhomogeneity-related horizontal transport (recalling low shear means less shear-production of TKE).
Beyond the comparison of derived analytical forms with measurements of conditional mean veer ⟨∆φ/∆z|α⟩, some general trends were also noted.For a given α ≳ 0.2, ⟨∆φ/∆z|α⟩ was larger offshore than for the onshore cases (though we remind that larger α are relatively rarer offshore compared to onshore conditions); this larger mean veer for a given α is due to the ABL depth h generally being lower offshore (see e.g.Liu and Liang, 2010, for offshore and onshore h).Perhaps counterintuitively, over the forested site the mean veer ⟨∆φ/∆z|α⟩ was smaller than other sites.As for the mean veer, for α ≳ 0.3 the long-term variability σ ∆φ/∆z|α was also found to be larger offshore; this may have impact on yaw error statistics, and may be the subject of future research.Analogous to σ α|U found in Kelly et al. (2014a) for shear, an empirical expression for the standard deviation of veer conditioned on wind speed (σ ∆φ/∆z|U ) was also found, with an approximately exponential decrease with speed.

Ongoing and future work
While the current work provided both theoretical meteorological relations and practical forms for veer in terms of shear, it did so without explicit treatment of buoyancy nor turbulent transport.Some relations including stabiilty within |S|/|G| in the shear contribution to veer were developed and tested; however these were not included here, as they did not offer improvement, are seen to be beyond the scope of the current work, and might also require stability effects to be explicitly incorporated within the cross-stress terms.Ongoing work involves addressing the latter: i.e., self-consistent α-based description of stability within the veer formulations, within both the shear and cross-stress contribtions in concert with the stability-perturbed geostrophic https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.drag law (Arya and Wyngaard, 1975;Kelly and Troen, 2016).Future work includes incorporation terrain-induced turbulent transport parameterization (following e.g.Kelly, 2020) into the veer, as well as study of the latter via LES.
Because the veer at commonly-occuring speeds (which occur below typical rated power) and also the mean veer are larger than for commonly-assumed neutral conditions, and since we have found relations for veer variability in terms of wind speed, practical ongoing work also involves vertical extrapolation of veer and accounting for its effect on power production.Accompanying this is validation and uncertainty quantification, towards pre-construction resource assessment as well as loads calculations.
https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.equation, which in flow-following coordinates at height z is solution, though one can note limits of the veer by considering two canonical cases where it can be solved: the Ekman and Ellison regimes, corresponding to simple prescriptions for ν T .Such limits were considered by van der Laan et al. (2020) for the geostrophic drag coefficient c G ≡ u * /|G| and ABL-integrated veer (cross-isobar angle) γ 0 ≡ φ 0 − φ G .

Figure 1 .
Figure 1.Veer behavior (plotted as degrees clockwise) for analytical/limiting cases of Ekman and Ellison.Left: veer versus shear exponent, for any Ekman or Ellison ABL depth h; Ekman is dotted purple, Ellison is dashed pink.Center/right: profiles of bulk veer for different ∆z;Ellison solution in right-hand plot is numerical solution (without approximation).
implemented into the k-ε turbulence closure equations solved by Ellipsys1D, as outlined in van der Laan et al. (2020); this includes use of small ambient values of turbulence intensity and dissipation rate above the ABL, with k-ε constants C µ = 0.03, C ε1 = 1.21,C ε2 = 1.92, σ k = 1.0 and σ ε = 1.3.The k-ε model provides the stresses occuring in the RANS equations, via the flux-gradient relation and ν T = C µ k 2 /ε; thus we see that such turbulence closure gives stresses aligned with velocity gradients.
https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.z 0 and ℓ max ; for lengthscale-limited k-ε RANS in the neutral ABL, one further has the relation h = ℓ 0the two Rossby numbers Ro 0 and Ro h for describing flow cases (van der Laan et al., 2020).Simulations were done over the full range of ABL depths, surface roughnesses, and wind speeds encountered in nature, which correspond to a range of Rossby numbers spanning 10 5 <Ro 0 < 10 10 and 15.8 <Ro h < 661.For simplicity |G| was set to 10 m s −1 and f to 10 −4 s −1 in the simulation set spanning these ranges of Rossby numbers.However, we remind that Rossby similarity means that for a given pair of {Ro 0 ,Ro h } and {z 0 , h} one has many (infinite) combinations of {|G|, f } which give the same |G|/f and thus the same dimensionless profile shapes of velocity, i.e. speed and direction as a function of dimensionless height zf /|G|.At any rate, the simulations cover ranges of (exceeding): ABL depths of 200-2000 m; roughness lengths from water's roughness (0.1 mm) up to 2.5 m; and |G| from 5-50 m s −1 .

Figure 3 .
Figure 3. Roughness dependence of turning (in degrees, clockwise) seen for two representative ABL depths, from 1D RANS simulations over a range of roughness lengths plotted directly against z0 (left) and alternately versus 1/ ln(Ro0) (right).Lines in right-hand plot indicate linear trend.
Figure 4 shows veer across three different rotor extents (z =50-100 m, 50-150 m, and 100-200 m), over a wide range of effective ABL depth h and associated Rossby number Ro h , for a commonly found roughness over land (1.6 cm).In Fig. 4 results are shown

Figure 4 .
Figure 4. Influence of ABL depth and associated Rossby number on veer (clockwise) for different turbine rotor spans.

Figure 5 .
Figure 5. Profiles of contributions to d cos γ/dz in (12) due to shear (left), stress gradients with Coriolis (center), and their sum (right).Five RANS simulations shown (two roughnesses and two ABL depths over land, one over sea) over typical turbine rotor heights; the listed z0 and h correspond to Rossby numbers using G = 10 m s −1 and f = 10 −4 s −1 .
Figure 6.Turning (bulk veer) in degrees clockwise versus shear exponent calculated from 50-150 m, over ranges of ABL depth and surface roughness; each point represents one RANS solution.Left: using G = 10 m s −1 and f = 10 −4 s −1 , over range of ABL depths spanning the values used in Figs.2-5 for water and typical land roughness.Center: again with G = 10 m s −1 and f = 10 −4 s −1 , over range of z0 spanning those used in Figs.2-5, for two ABL (typical) depths used in previous figures.Right: over wider range of {Ro0, Ro h } spanning that found in nature; note larger vertical axis scale.
https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.1997), from 80-200 m height for two years; and one year from a commercial site dubbed 'MR' which sits on a ridge over a mostly forested (>∼ 3/4) area but dominated by hills having elevation differences up to ∼200 m within 10 km distance, using anemometers at 40-136 m height. 11

Figure 8 .
Figure 8. Distributions of shear exponent (left) and corresponding veer (right) at Høvsore between 60-160 m, from the homogeneous eastern land sectors (red) and from the sea sectors to the west (blue).Solid lines are for measurements spanning 60-160 m; dashed for those spanning 100-160 m.
vary with height more for the offshore flow than for the homogeneous land directions.The change of ⟨α⟩ from 60-160 m to https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.

Figure 11 .
Figure11.Joint distribution of veer and wind speed at sites considered.Solid line shows ⟨∆φ/∆z⟩ U , calculated using 1 m s −1 bins; lightest shades are 2% as likely as darkest color in each plot.

Figure 12 .
Figure 12.Measured standard deviation of veer conditioned on wind speed, again using 1 m s −1 bins, for the sites considered.Neutral conditions at Høvsøre defined as in earlier figures, i.e. |L −1 | < 0.001 m −1 .

Figure 13 .
Figure 13.Top: joint distribution of veer and shear exponent observed over 10 years from 60-160 m for the Høvsøre land sectors in neutral conditions, in different speed ranges; axes zoomed in to show detail, and occurence rate normalized per wind speed range (each plot has a different color scale, showing occurence rate in increments of 1/10, with lightest representing 10% as likely as darkest shade).Bottom: the same joint distributions shown with unscaled rate of occurrence (number of counts per {α, ∆φ/∆z} bin); axis ranges are chosen to compare with later figures.

Figure 14 .
Figure 14.Joint distribution of veer and shear exponent for different speed ranges from Høvsøre land sectors, but for all stability conditions; plots are analogous to those in bottom of Fig. 13.All plots use same color scale; color bar denotes count.
Fig. 15), while for the more stability-dominated homogeneous cases (top plots in Fig. 15) a value of c sα =0.7 for Høvsøre and 0.8 for 595 Cabauw gave reasonable fits.The latter aspect could be practically addressed by directly casting c sα as a minimal value plus an amount depending on the long-term variability in positive stability (labelled σ + following Kelly and Gryning, 2010); we https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.note Cabauw has larger values of σ + than Høvsøre, which has larger σ + then Østerild.However, obtaining such an expression is beyond the scope of the current article, and some sites could have factors other than stability which enhance the veer.We do find that including stability within the drag law via M-O theory (for positive L −1 values consistent with observed distributions)

Figure 16 .
Figure 16.Statistics of veer conditioned on shear exponent, across sites considered.
https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.roughness z 0 and ABL depth h is also elucidated (through Rossby numbers Ro 0 and Ro h defined by each), through numerical solution of the 1-D RANS equations under neutral conditions with lengthscale-limited k-ε turbulence closure (i.e.neutral but also translatable to stable conditions, see van der Laan et al., 2020).
https://doi.org/10.5194/wes-2022-119Preprint.Discussion started: 9 January 2023 c Author(s) 2023.CC BY 4.0 License.with the capping inversion, ∆φ/∆z increases overall with z through the jet as the surface-based stress decreases with height (though there can be occasional deviations from this behavior due to stress profiles affected by upwind inhomogeneities or large coherent structures).