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 for vertical extrapolation of mean wind speeds, as well as being a key parameter for wind turbine load 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 crosswind stress. Following from the theoretical derivations, which are neither limited to the surface layer nor 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 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 in veer conditioned on wind speed. Analyzing observed joint probability distributions of α 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.e., conditioned on) shear exponent; lastly, the limitations, applicability, and behavior of these forms are discussed along with their use and further developments for both meteorology and wind energy.

1 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 this (Triviño, 2017; Kelly et al., 2019b). Shear is also a key parameter for flow characterization towards load 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 labeled 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 they 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 6-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.5 MW) clockwise-turning turbine rotor to reduce power production, opposite to 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), who 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 effects (Hulsman et al., 2022; Narasimhan et al., 2022).

In this paper we investigate wind veer, showing its joint behaviors with and connection to shear and key parameters used to describe atmospheric boundary layer flow. In Sect. 2, after reviewing the expression of the shear exponent and its relation to stability and turbulence, we derive new relations for veer; we show veer to be composed of shear-driven and Coriolis-associated stress gradient contributions. The theoretical behavior of veer is also derived for canonical cases such as Ekman and surface-layer flow, as well as the effect of shear-stress misalignment on veer. Further, in Sect. 2.4 practical relations from micrometeorology are elucidated towards an evaluation of the expressions developed for veer. Section 3 includes an analysis of veer, exploring and connecting the developed relations to both computational modeling and observations. Section 3.1 gives Reynolds-averaged Navier–Stokes (RANS; mean) simulation results over flat terrain in neutral conditions for hundreds of combinations of surface Rossby number and ABL-depth Rossby number, showing the dependence of veer on the latter, as well as the counteracting behavior of veer's two primary components. Section 3.2 begins with an analysis of multi-year observations from six different flow regimes across four sites showing the statistical behavior of shear with stability and subsequently that of veer, also providing new empirical relations for the probability of occurrence of larger veer (due to the effect of stable conditions) and for the variability in veer with wind speed. The observational analysis concludes in Sect. 3.3 with simplified practical relations for veer based on observed shear, including a comparison with joint distributions of veer and shear across the six flows analyzed. Finally the results are summarily discussed and conclusions given, with ongoing and future work also described for the reader.

2 Theory and development

In this section we define the shear exponent and veer and then derive relations for veer in terms of shear and vertical gradients of stress, as mentioned in the previous paragraph. Section 2.3 provides a number of expressions for veer; this is done to facilitate its calculation and interpretation in the different coordinate systems typically considered in wind energy flow analyses, and we also include forms that are independent of coordinate system. Because coordinates aligned with the mean wind for a given height of interest (e.g., hub height) are commonly used in wind energy and because expressions for veer in such a coordinate system are simpler to express and calculate, we ultimately arrive at two forms in such a system (Eqs. 14 and 16); due to its robustness, one of these (Eq. 14) will later be shown in Sect. 3.3 to be further simplifiable and usable (as Eqs. 39 or 40) in comparison with measurements.

2.1 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,

$$\begin{array}{}\text{(1)}& \mathit{\alpha }\equiv {\displaystyle \frac{\partial U/\partial z}{U/z}}={\displaystyle \frac{\partial \ln U}{\partial \ln z}},\end{array}$$

and the dimensionless wind speed gradient,

$$\begin{array}{}\text{(2)}& {\mathrm{\Phi }}_m\equiv {\displaystyle \frac{\mathrm{d}U/\mathrm{d}z}{u_{*\mathrm{0}}/\mathit{\kappa }z}}={\displaystyle \frac{\mathit{\kappa }U}{u_{*\mathrm{0}}}}\mathit{\alpha },\end{array}$$

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¹. Note that Eq. (1) is derived from the power-law expression for wind speed,

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{U}{U_{\mathrm{ref}}}}={\left({\displaystyle \frac{z}{z_{\mathrm{ref}}}}\right)}^{\mathit{\alpha }},\end{array}$$

which is assumed to be valid over some extent around height z_ref, with U_ref≡U(z_ref). The power-law (Eq. 3) with shear exponent (Eq. 1) has been used in wind engineering for decades (e.g., Irwin, 1979; Petersen et al., 1998) due to its simplicity and because it does not require any information other than the wind speed at two heights. Although Eqs. (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–Obukhov (M–O) theory Φ_m is a function of the stability $z/L$ which is proportional to surface heat flux H₀ divided by $u_{*\mathrm{0}}^{\mathrm{3}}$, i.e., the reciprocal Obukhov length is $\mathrm{1}/L=\mathit{\kappa }(g/T_{\mathrm{0}})H_{\mathrm{0}}/u_{*\mathrm{0}}^{\mathrm{3}}$, where T₀ is the background temperature and g is the gravitational acceleration (Monin and Obukhov, 1954); the Φ_m function and corresponding M–O wind profile (which arises via integrating $\mathrm{d}U/\mathrm{d}z$ in Eq. (2) from a height equal to the roughness length z₀ up to height z) thus require a number of assumptions and more information than the calculation of α via Eq. (1) or use of the power-law (Eq. 3). Monin–Obukhov wind profiles also require the surface roughness length (z₀), while the friction velocity u_*0 (and thus shear stress) is assumed to be constant in the surface layer where M–O theory is most valid²; further, the assumptions of stationarity and a uniform flat surface are implicit in the use of M–O theory. Following surface-layer theory one could write an equivalent shear exponent ${\mathit{\alpha }}_{\mathrm{ASL}}={\mathrm{\Phi }}_m(z/L)/\left[\ln \right(z/z_{\mathrm{0}})-{\mathrm{\Psi }}_m(z/L\left)\right]$, where

$${\mathrm{\Psi }}_m=\underset{z_{\mathrm{0}}}{\overset{z}{\int }}[\mathrm{1}-{\mathrm{\Phi }}_m(z^{\prime }/L\left)\right]\mathrm{d}\ln z^{\prime }$$

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³.

2.1.1 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

$$\begin{array}{}\text{(4)}& \mathit{\alpha }={\displaystyle \frac{z}{U}}{\displaystyle \frac{(\mathit{\epsilon }-B-T)}{-\langle uw\rangle }}\end{array}$$

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. We point out that Kelly et al. (2014a) ignored crosswind stress 〈vw〉 when deriving Eq. (4); however, it still shows that, e.g., shear will increase in stable conditions (B<0) and decrease in unstable conditions (B>0), as will be demonstrated using observations in Sect. 3.2; further, as we will see in Sect. 2.3, this is also related to the veer. Within the ASL under these conditions where M–O theory is valid and 〈vw〉→0, using the neutral value of dissipation rate as ${\mathit{\epsilon }}_{\mathrm{0}}\equiv u_{*\mathrm{0}}^{\mathrm{3}}/\left(\mathit{\kappa }z\right)$ along with the dimensionless functions ${\mathrm{\Phi }}_{\mathit{\epsilon }}\equiv \mathit{\epsilon }/{\mathit{\epsilon }}_{\mathrm{0}}$ and ${\mathrm{\Phi }}_T\equiv T/{\mathit{\epsilon }}_{\mathrm{0}}$ (Kaimal and Finnigan, 1994), we can express a surface-layer version of Eq. (4) as

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\mathrm{ASL}}={\displaystyle \frac{u_{*\mathrm{0}}}{\mathit{\kappa }U}}\left({\mathrm{\Phi }}_{\mathit{\epsilon }}+{\displaystyle \frac{z}{L}}-{\mathrm{\Phi }}_T\right)\approx I_u\left({\mathrm{\Phi }}_{\mathit{\epsilon }}+{\displaystyle \frac{z}{L}}-{\mathrm{\Phi }}_T\right),\end{array}$$

since by definition $B/{\mathit{\epsilon }}_{\mathrm{0}}=-z/L$ and $u_{*\mathrm{0}}^{\mathrm{2}}=-\langle uw\rangle$; here $I_u\equiv {\mathit{\sigma }}_u/U$ is the streamwise turbulence intensity. The dimensionless dissipation rate (M–O function) Φ_ε≥1 is roughly $\mathrm{1}+\mathrm{5}z/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^{-\mathrm{1}}>\mathrm{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 Eq. (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.

2.2 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\equiv U+iV=\left|S\right|e^{i\mathit{\phi }}$. For a mean wind direction defined at some height z, the veer can be defined as a directional shear $\partial \mathit{\phi }/\partial z$ through the wind direction:

$$\begin{array}{}\text{(6)}& \mathit{\phi }\left(z\right)=\arg \left[S\right(z\left)\right]=\arctan \left[{\displaystyle \frac{V\left(z\right)}{U\left(z\right)}}\right].\end{array}$$

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 ${\mathit{\phi }}_{\mathrm{0}}\equiv \mathit{\phi }\left(\mathrm{0}\right)=\mathrm{0}$). We follow this convention unless stated otherwise, as done for some expressions later in Sect. 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 $\left|G\right|$, and surface roughness length z₀ 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 $\mathrm{\Delta }z/h$ of the ABL will also depend on these parameters; thus for a given site and height, $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial S}{\partial t}}=\mathrm{0}=-if(S-G)-{\displaystyle \frac{\partial \langle sw\rangle }{\partial z}}\end{array}$$

for stationary and horizontally homogeneous conditions (thus neglecting advection). Here the kinematic horizontal pressure gradient $\mathrm{\nabla }p/\mathit{\rho }=f\mathit{\{}{V}_{\mathrm{G}},-U_{\mathrm{G}}\mathit{\}}$ is also written like a velocity in complex form as $G\equiv U_{\mathrm{G}}+i{V}_{\mathrm{G}}=(-\partial p/\partial y+i\partial p/\partial x)/\left(\mathit{\rho }f\right)$. The mean stresses are dominated by vertical momentum transport 〈sw〉, where w denotes (turbulent) vertical velocity fluctuations and $s\equiv u+iv$ the horizontal velocity fluctuations.

At a given height z, taking the differential of Eq. (6) (recalling $\mathrm{d}\arctan x=\mathrm{d}x/[\mathrm{1}+x^{\mathrm{2}}]$ and using the chain rule) gives

$$\begin{array}{}\text{(8)}& \mathrm{d}\mathit{\phi }={\displaystyle \frac{U\mathrm{d}V-V\mathrm{d}U}{|S{|}^{\mathrm{2}}}};\end{array}$$

here the superscript asterisk denotes a complex conjugate. Applying $\partial /\partial z$ to Eqs. (7) and (8) and combining provides a basic expression for veer:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}={\displaystyle \frac{U}{|S{|}^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{1}}{f}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\langle uw\rangle }{\partial z^{\mathrm{2}}}}+{\displaystyle \frac{\partial V_{\mathrm{G}}}{\partial z}}\right]+{\displaystyle \frac{V}{|S{|}^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{1}}{f}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\langle vw\rangle }{\partial z^{\mathrm{2}}}}-{\displaystyle \frac{\partial U_{\mathrm{G}}}{\partial z}}\right].\end{array}$$

In the case of zero geostrophic shear ($\mathrm{d}\mathbf{G}/\mathrm{d}z=\mathrm{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 Eq. (9) can be written more simply as

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\phi }\left(z\right)}{\partial z}}{|}_{\mathrm{d}\mathbf{G}/\mathrm{d}z\to \mathrm{0}}={\displaystyle \frac{\mathrm{1}}{f\left|S\right|}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\langle uw\rangle {|}_{{\mathbf{e}}_x\parallel \mathbf{U}\left(z\right)}.\end{array}$$

Though Eqs. (9) and (10) are not directly very useful for relating veer to shear, since the shear is implicit in the stress terms (and one would need to know the profiles of horizontal stresses to use these equations), they do illustrate that the curvature of stress profiles and Coriolis effect are the basis for mean veer following Eq. (7) and also that geostrophic shear can further contribute to veer (e.g., due to baroclinity; see Hoxit, 1974; Arya and Wyngaard, 1975; Pedersen et al., 2013).

2.3 Relating veer to shear

Towards relating the veer to shear, one can alternately derive the veer by first taking the time derivative of Eq. (8); using the real and imaginary parts of Eq. (7), in the horizontally homogeneous limit (ignoring advection) one obtains a rate equation for mean wind direction:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial t}}=\left[{\displaystyle \frac{V}{|S{|}^{\mathrm{2}}}}{\displaystyle \frac{\partial \langle uw\rangle }{\partial z}}-{\displaystyle \frac{U}{|S{|}^{\mathrm{2}}}}{\displaystyle \frac{\partial \langle vw\rangle }{\partial z}}\right]+f\left({\displaystyle \frac{\left|G\right|}{\left|S\right|}}\cos \mathit{\gamma }-\mathrm{1}\right).\end{array}$$

The “turning” angle $\mathit{\gamma }\equiv \mathit{\phi }-{\mathit{\phi }}_{\mathrm{G}}$ between geostrophic and mean wind directions (e.g., Wyngaard, 2010) arises through⁴

$$U_{\mathrm{G}}U+V_{\mathrm{G}}V=\mathbf{U}\cdot \mathbf{G}=\left|S\right|\left|G\right|\cos \mathit{\gamma },$$

by taking $\partial /\partial t$ of Eq. (6) or equivalently $U\partial V/\partial t-V\partial U/\partial t$ via Eq. (8). The geostrophic wind direction is defined as ${\mathit{\phi }}_{\mathrm{G}}\equiv \arctan (V_{\mathrm{G}}/U_{\mathrm{G}})$, and the “cross-isobar” angle, i.e., the turning over the whole ABL $({\mathit{\gamma }}_{\mathrm{0}}={\mathit{\phi }}_{\mathrm{0}}-{\mathit{\phi }}_{\mathrm{G}})$, is generally less than 45^∘ (Grisogono, 2011)⁵; 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⁶. Note that φ, and thus γ, can vary with height z (as can φ_G in baroclinic conditions).

Assuming statistical stationarity so that $\partial \mathit{\phi }/\partial t=\mathrm{0}$, the vertical derivative of Eq. (11) can be written most conveniently in terms of the deviation of dimensionless wind from streamwise; taking the vertical derivative of Eq. (11) if we again take $\mathrm{d}\mathbf{G}/\mathrm{d}z=\mathrm{0}$ (neglect baroclinity), then

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial \cos \mathit{\gamma }}{\partial z}}={\displaystyle \frac{\mathrm{1}}{\left|G\right|}}{\displaystyle \frac{\partial \left|S\right|}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{f\left|G\right|}}{\displaystyle \frac{\partial }{\partial z}}\left[{\displaystyle \frac{U}{\left|S\right|}}{\displaystyle \frac{\partial \langle vw\rangle }{\partial z}}-{\displaystyle \frac{V}{\left|S\right|}}{\displaystyle \frac{\partial \langle uw\rangle }{\partial z}}\right].\end{array}$$

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. Equation (12) clearly separates the shear and stress–Coriolis 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, from hub to tip, or between measurement and hub heights. Re-expressing Eq. (11) with the coordinate system defined by having x in the mean wind direction at height z so that S(z)=U(z)e_x and $\left|S\right(z\left)\right|=U\left(z\right)$, in the mean (for $\mathrm{d}\mathit{\phi }/\mathrm{d}t=\mathrm{0}$) one has

$$\begin{array}{}\text{(13)}& \cos \mathit{\gamma }={\displaystyle \frac{\left|S\right|}{\left|G\right|}}+{\displaystyle \frac{\mathrm{1}}{f\left|G\right|}}{\displaystyle \frac{\partial \langle vw{\rangle }_{⟂}}{\partial z}},\end{array}$$

where we use the shorthand notation 〈vw〉_⟂ to denote the stress perpendicular to the mean flow at a given height. Taking the inverse cosine and subsequently the vertical derivative, noting that $\partial \mathit{\gamma }/\partial z=\partial \mathit{\phi }/\partial z$ and $\mathrm{d}\left(\arccos x\right)=-\mathrm{d}x/\sqrt{\mathrm{1}-x^{\mathrm{2}}}$ while recalling $\partial \left|S\right|/\partial z=\mathit{\alpha }\left|S\right|/z$, we get

$$\begin{array}{ll}{\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{|}_{{\mathbf{e}}_x\parallel \mathbf{U}\left(z\right)}=& {\displaystyle }-\left[{\displaystyle \frac{\left|S\right|}{\left|G\right|}}{\displaystyle \frac{\mathit{\alpha }}{z}}+{\displaystyle \frac{\mathrm{1}}{f\left|G\right|}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\langle vw{\rangle }_{⟂}}{\partial z^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\left|S\right|}}{\displaystyle \frac{\partial V}{\partial z}}{\displaystyle \frac{\partial \langle uw{\rangle }_{\parallel }}{\partial z}}\right)\right]\\ {\displaystyle }& {\displaystyle }{\left[\mathrm{1}-{\left({\displaystyle \frac{\left|S\right|}{\left|G\right|}}+{\displaystyle \frac{\mathrm{1}}{f\left|G\right|}}{\displaystyle \frac{\partial \langle vw{\rangle }_{⟂}}{\partial z}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}};\end{array}$$

here the subscript ∥ is used to remind us that 〈uw〉 is parallel to the flow at height z. Further, this coordinate system (Eq. 8) gives $\partial V/\partial z=U\partial \mathit{\phi }/\partial z=\left|S\right|\partial \mathit{\phi }/\partial z$ so that the $\partial V/\partial z$ part just becomes an additional veer term on the right-hand side; collecting the $\partial \mathit{\phi }/\partial z$ on the left side and rearranging we then obtain

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{|}_{{\mathbf{e}}_x\parallel \mathbf{U}\left(z\right)}={\displaystyle \frac{-\frac{\left|S\right|}{\left|G\right|}\frac{\mathit{\alpha }}{z}-\frac{\mathrm{1}}{f\left|G\right|}\frac{{\partial }^{\mathrm{2}}\langle vw{\rangle }_{⟂}}{\partial z^{\mathrm{2}}}}{{\left[\mathrm{1}-{\left(\frac{\left|S\right|}{\left|G\right|}+\frac{\mathrm{1}}{f\left|G\right|}\frac{\partial \langle vw{\rangle }_{⟂}}{\partial z}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}-\frac{\mathrm{1}}{f\left|G\right|}\frac{\partial \langle uw{\rangle }_{\parallel }}{\partial z}}}.\end{array}$$

As will be demonstrated in a later section, basically one sees from the numerator of Eq. (14) that the veer is comprised of a shear-associated part and a crosswind stress-curvature part; the denominator is basically 1 minus a few relatively small terms. The more generic form of veer, for an arbitrary coordinate system, also follows from Eq. (11):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}={\displaystyle \frac{-\frac{\left|S\right|}{\left|G\right|}\frac{\mathit{\alpha }}{z}+\frac{\mathrm{1}}{f\left|G\right|}\frac{\partial }{\partial z}\left(\frac{V}{\left|S\right|}\frac{\partial \langle uw\rangle }{\partial z}-\frac{U}{\left|S\right|}\frac{\partial \langle vw\rangle }{\partial z}\right)}{\sqrt{\mathrm{1}-{\left[\frac{\left|S\right|}{\left|G\right|}-\frac{\mathrm{1}}{f\left|G\right|}\left(\frac{V}{\left|S\right|}\frac{\partial \langle uw\rangle }{\partial z}-\frac{U}{\left|S\right|}\frac{\partial \langle v\rangle }{\partial z}\right)\right]}^{\mathrm{2}}}}}.\end{array}$$

We note that Eqs. (14) and (15) avoid the use of the turning (ageostrophic) angle γ, and subsequently nonlinear functions involving φ_G, which becomes apparent if one expands cos γ in Eqs. (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 Eqs. (12) and (15), $U/\left|S\right|$ and $V/\left|S\right|$ can be written as cos φ and sin φ, respectively. Then from Eq. (12) and using $\cos \mathit{\gamma }=\cos \mathit{\phi }\cos {\mathit{\phi }}_{\mathrm{G}}+\sin \mathit{\phi }\sin {\mathit{\phi }}_{\mathrm{G}}$, again in coordinates defined by $\left|S\right(z\left)\right|=U\left(z\right)$, after some rearranging we arrive at an expression for veer like Eq. (14):

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{|}_{{\mathbf{e}}_x\parallel \mathbf{U}\left(z\right)}={\displaystyle \frac{\frac{\left|S\right|\mathit{\alpha }}{\left|G\right|z}+\frac{{\partial }^{\mathrm{2}}\langle vw{\rangle }_{⟂}/\partial z^{\mathrm{2}}}{f\left|G\right|}}{\sin {\mathit{\phi }}_{\mathrm{G}}+\frac{\partial \langle uw{\rangle }_{\parallel }/\partial z}{f\left|G\right|}}}.\end{array}$$

Compared to Eq. (14) this lacks a negative sign, but sin φ_G is negative and with a larger magnitude than the positive contribution to the denominator, $\partial \langle uw{\rangle }_{\parallel }/\partial z/\left(f\right|G\left|\right)$; this will become more apparent in the sections which follow. We also note that in these coordinates φ_G=γ_G(z), and opposite signs will occur for the Southern Hemisphere (Eqs. 14–16 give $\mathrm{d}\mathit{\phi }/\mathrm{d}z$ signed for the Northern Hemisphere in mathematical coordinates, i.e., negative), reflecting winds rotating on average clockwise with increasing height.

For wind energy $\partial \left(\cos \mathit{\phi }\right)/\partial z$ might be considered relevant as $\partial \mathit{\phi }/\partial z$ because it allows the direct expression of the veer-induced variation in the streamwise wind velocity component relative to a reference height such as hub height. One could expect that the reduction in 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 Eq. (12).

One last relation between shear and veer can also be elucidated by considering a corrected version of Eq. (4). By keeping the lateral shear term $\langle vw\rangle \partial V/\partial z$ in the TKE rate equation and then again using coordinates defined with x in the mean direction at height z and subsequently $\partial V/\partial z\to U\partial \mathit{\phi }/\partial z$, then Eq. (4) contains an additional contribution, becoming

$$\begin{array}{}\text{(17)}& \mathit{\alpha }{|}_{{\mathbf{e}}_x\parallel \mathbf{U}\left(z\right)}={\displaystyle \frac{\mathit{\epsilon }-B-T}{-U\langle uw{\rangle }_{\parallel }/z}}-z{\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}{\displaystyle \frac{\langle vw{\rangle }_{⟂}}{\langle uw{\rangle }_{\parallel }}}.\end{array}$$

Recalling in the ABL that $\langle uw{\rangle }_{\parallel }<\mathrm{0}$ (momentum gets transferred towards the surface⁷) because $\langle vw{\rangle }_{⟂}>\mathrm{0}$ in the ABL (Wyngaard, 2010), we see as in Eqs. (14)–(16) that negative $\partial \mathit{\phi }/\partial z$ (clockwise veer) is associated with positive shear; note that the sign of $\partial \mathit{\phi }/\partial z$ is flipped in typical wind energy coordinates (left-handed, with 0^∘ corresponding to wind from the north and increasing clockwise). Although we have provided Eq. (17) to both improve Eq. (4) from Kelly et al. (2014a) and offer insight into how shear and veer are linked within the context of TKE, we advise that it is not easily utilized compared to forms like Eq. (14); the latter will be applied and investigated further in later sections.

2.3.1 Misalignment of shear and stress

One can see a connection between the shear, veer, and stress in Eqs. (9) and (12), and we can further examine the relation between shear and stress using complex notation as in Eq. (7). The “misalignment” can be expressed via the angle between $\partial S/\partial z$ and 〈sw〉, i.e.,

$$\begin{array}{}\text{(18)}& {\mathit{\beta }}_{\mathrm{ma}}\equiv \left(\mathit{\phi }-{\mathit{\phi }}_{sw}\right)=\arg (\partial S/\partial z)-\arg (\langle sw\rangle ).\end{array}$$

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 $\partial \langle uw\rangle /\partial t$ to $i\partial \langle vw\rangle /\partial t$, we may write

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial \langle sw\rangle }{\partial t}}=\mathrm{0}\simeq \langle w^{\mathrm{2}}\rangle {\displaystyle \frac{\partial S}{\partial z}}-{\displaystyle \frac{\langle sw\rangle }{{\mathit{\tau }}_{\mathrm{R}}}}-{\displaystyle \frac{\partial }{\partial z}}\langle sww\rangle .\end{array}$$

The pressure-strain contribution has been written as $\langle sw\rangle /{\mathit{\tau }}_{\mathrm{R}}$ via the commonly used Rotta (1951) parameterization, where τ_R is the Rotta timescale; this is the basis for commonly used flux-gradient relations (Wyngaard, 2004). In such mixing-length relations, i.e., using the “Boussinesq hypothesis”, 〈w²〉τ_R is simply written as a turbulent diffusivity −ν_T, and the final term in Eq. (19) is neglected. We continue to neglect advection and horizontal transport (such as $U\partial \langle sw\rangle /\partial x$ and $\partial \langle suw\rangle /\partial x$, respectively); these can also contribute to misalignment between $\partial S/\partial 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 two-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 provide a convenient description of the veer, it can be re-cast to show the effect of the misalignment angle:

$$\begin{array}{}\text{(20)}& f\left|G\right|\sin \mathit{\gamma }=-{\displaystyle \frac{\partial |\langle sw\rangle |}{\partial z}}\cos {\mathit{\beta }}_{\mathrm{ma}}-|\langle sw\rangle |\sin {\mathit{\beta }}_{\mathrm{ma}}{\displaystyle \frac{\partial {\mathit{\phi }}_{sw}}{\partial z}}.\end{array}$$

Thus when the stress is aligned with the shear (β_ma=0), then $f\left|G\right|\sin \mathit{\gamma }=-\partial |\langle sw\rangle |/\partial z$; this can be seen as a case of Eq. (13). The contribution of stress–shear misalignment to the veer can also be seen considering Eq. (19) with our earlier derivations, with misalignment modifying the stresses. For example the crosswind stress in Eqs. (13)–(15) can be written

$$\begin{array}{}\text{(21)}& \langle vw{\rangle }_{⟂}=-{\mathit{\nu }}_T\left[{\displaystyle \frac{\partial V_{⟂}}{\partial z}}+{\displaystyle \frac{\partial \langle vww{\rangle }_{⟂}/\partial z}{\mathrm{2}k/\mathrm{3}}}\right],\end{array}$$

since the Rotta timescale can be expressed in terms of turbulent kinetic energy k via ${\mathit{\nu }}_{\mathrm{T}}={\mathit{\tau }}_{\mathrm{R}}\langle uu+vv+ww\rangle /\mathrm{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.

2.3.2 Canonical solutions using an eddy diffusivity

When turbulent transport of stress is negligible (along with baroclinity and inhomogeneity), in steady conditions the stress and mean velocity gradient are aligned. This allows the use of an eddy diffusivity ν_T to express the stress as $\langle sw\rangle =-{\mathit{\nu }}_{\mathrm{T}}\left(z\right)\partial S/\partial z$ which can then be cast as a nonlinear differential equation using the stress cast in terms of eddy diffusivity and shear in Eq. (10), which again in flow-following coordinates at height z neglecting geostrophic shear is

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\partial {\mathit{\phi }}_{{\mathit{\nu }}_{\mathrm{T}}}}{\partial z}}={\displaystyle \frac{-\mathrm{1}}{f\left|G\right|}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\left[{\mathit{\nu }}_{\mathrm{T}}\left(z\right){\displaystyle \frac{\partial U}{\partial z}}\right].\end{array}$$

This defies the solution without some prescription for ν_T(z), 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_{\mathrm{G}}\equiv u_{*}/\left|G\right|$ and ABL-integrated veer (cross-isobar angle) ${\mathit{\gamma }}_{\mathrm{0}}\equiv {\mathit{\phi }}_{\mathrm{0}}-{\mathit{\phi }}_{\mathrm{G}}$.

2.3.3 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 $\langle sw\rangle =-{\mathit{\nu }}_{\mathrm{Ek}}\partial S/\partial z$. Thus the momentum balance (Eq. 7) simplifies to

$$\begin{array}{}\text{(23)}& f\left(S_{\mathrm{Ek}}-G\right)=-i{\mathit{\nu }}_{\mathrm{Ek}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{S}_{\mathrm{Ek}}}{\partial z^{\mathrm{2}}}},\end{array}$$

which gives the classic Ekman solution:

$$\begin{array}{}\text{(24)}& S_{\mathrm{Ek}}=G\left(\mathrm{1}-e^{-(\mathrm{1}+i)z/h_{\mathrm{Ek}}}\right),\end{array}$$

where the characteristic Ekman (e-folding) height h_Ek is defined as $h_{\mathrm{Ek}}\equiv \sqrt{\mathrm{2}{\mathit{\nu }}_{\mathrm{Ek}}/f}$. Simpler than relating Ekman veer to shear, the solutions above along with Eq. (9) give the veer directly (again in radians) as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\phi }}_{\mathrm{Ek}}}{\partial z}}& ={\displaystyle \frac{-{\mathit{\nu }}_{\mathrm{Ek}}}{f|S_{\mathrm{Ek}}{|}^{\mathrm{2}}}}\left[U{\displaystyle \frac{{\partial }^{\mathrm{3}}U}{\partial z^{\mathrm{3}}}}+V{\displaystyle \frac{{\partial }^{\mathrm{3}}V}{\partial z^{\mathrm{3}}}}\right]\\ & ={\displaystyle \frac{e^{-z/h_{\mathrm{Ek}}}}{h_{\mathrm{Ek}}}}\cdot {\displaystyle \frac{\cos \left(z/h_{\mathrm{Ek}}\right)-\sin \left(z/h_{\mathrm{Ek}}\right)-e^{-z/h_{\mathrm{Ek}}}}{\mathrm{1}-\mathrm{2}{e}^{-z/h_{\mathrm{Ek}}}\cos \left(z/h_{\mathrm{Ek}}\right)+e^{-\mathrm{2}z/h_{\mathrm{Ek}}}}}\\ & \simeq {\displaystyle \frac{-\mathrm{0.5}+z/\left(\mathrm{6}{h}_{\mathrm{Ek}}\right)}{h_{\mathrm{Ek}}}};\end{array}\end{array}$$

this result has units of radians per meter measured counter-clockwise, with the linear approximation⁸ deviating from the exact form by less than 1 % for z<1.5h_Ek. Integrated over $z\pm \mathrm{\Delta }z/\mathrm{2}$, this gives the veer across an extent Δz:

$$\begin{array}{}\text{(26)}& \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{Ek}}\simeq {\displaystyle \frac{-\mathrm{\Delta }z}{\mathrm{2}{h}_{\mathrm{Ek}}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }z}{\mathrm{3}{h}_{\mathrm{Ek}}}}\right).\end{array}$$

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 γ₀.

From Eq. (24) one can also find an expression for the Ekman shear exponent α_Ek via Eq. (1):

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{Ek}}& ={\displaystyle \frac{|\partial S_{\mathrm{Ek}}/\partial z|}{\left|S_{\mathrm{Ek}}\right|/z}}={\displaystyle \frac{\sqrt{\mathrm{2}}(z/h_{\mathrm{Ek}})}{\sqrt{\mathrm{1}-\mathrm{2}\cos \left(z/h_{\mathrm{Ek}}\right)e^{z/h_{\mathrm{Ek}}}+e^{\mathrm{2}z/h_{\mathrm{Ek}}}}}}\\ & \simeq \left(\mathrm{1}-{\displaystyle \frac{\sqrt{\mathrm{2}}}{\mathit{\pi }}}{\displaystyle \frac{z}{h_{\mathrm{Ek}}}}\right).\end{array}\end{array}$$

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.

Figure 1 Veer behavior (plotted as degrees clockwise) for analytical and limiting cases of Ekman and Ellison. (a) Veer versus shear exponent for any Ekman or Ellison ABL depth h; Ekman is dotted purple, and Ellison is dashed pink. (b, c) Profiles of bulk veer for different Δz; Ellison solution in (c) is the numerical solution (without approximation).

2.3.4 Ellison solution (linear diffusivity profile and surface-layer regime)

Using a surface-layer eddy-viscosity relation ${\mathit{\nu }}_{\mathrm{T}}\left(z\right)=\mathit{\kappa }{u}_{*}z$ consistent with ASL theory, Ellison (1956) derived the solution of Eq. (7) for the (complex) wind vector, resulting in a profile of geostrophic “deficit” expressible as (Krishna, 1980)

$$\begin{array}{}\text{(28)}& G-S_{\mathrm{Ell}}\left(z\right)={\displaystyle \frac{\mathrm{2}{u}_{*}}{\mathit{\kappa }}}\left[{\ker }_{\mathrm{0}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}fz}{\mathit{\kappa }{u}_{*}}}}\right)+i{\mathrm{kei}}_{\mathrm{0}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}fz}{\mathit{\kappa }{u}_{*}}}}\right)\right],\end{array}$$

where ker₀(x) and kei₀(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 Eq. (24) with a complex argument, as

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{S}_{\mathrm{Ell}}\left(z\right)& =G-{\displaystyle \frac{\mathrm{2}{u}_{*}}{\mathit{\kappa }}}{K}_{\mathrm{0}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}ifz}{\mathit{\kappa }{u}_{*}}}}\right)\\ & =G\left[\mathrm{1}-\mathrm{2}{\displaystyle \frac{c_{\mathrm{G}}}{\mathit{\kappa }}}{K}_{\mathrm{0}}\left({\displaystyle \frac{(\mathrm{1}+i)z}{\sqrt{{\mathit{\nu }}_{\mathrm{T}}\left(z\right)/f}}}\right)\right]\\ & =G\left[\mathrm{1}-\mathrm{2}{\displaystyle \frac{c_{\mathrm{G}}}{\mathit{\kappa }}}{K}_{\mathrm{0}}\left(\sqrt{{\displaystyle \frac{\mathrm{2}iz}{h_{\mathrm{mE}}}}}\right)\right];\end{array}\end{array}$$

K₀(x) is the zeroth-order modified Bessel function of the second kind, and the modified Ekman length scale is defined by $h_{\mathrm{mE}}\equiv \mathit{\kappa }{u}_{*}/f$, also equal to ${\mathit{\nu }}_{\mathrm{T}}\left(z\right)/fz$. For the range $\mathrm{0.02}\lesssim c_{\mathrm{G}}\lesssim \mathrm{0.06}$ encountered in nature under neutral conditions (Hess and Garratt, 2002; van der Laan et al., 2020), for $z_H/h_{\mathrm{mE}}\gg \mathrm{0.1}$ the arctangent of $\mathrm{\Im }\mathit{\left\{}S\mathit{\right\}}/\mathrm{\Re }\mathit{\left\{}S\mathit{\right\}}$ can be approximated via series expansions of Eqs. (28) or (29) to yield the practical result

$$\begin{array}{}\text{(30)}& \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{Ell}}\left(z_H,\mathrm{\Delta }z\right)\approx \mathit{\pi }{c}_{\mathrm{G}}\exp \left(-\sqrt{z^{\prime }/z_H}\right){|}_{z_H-\mathrm{\Delta }z/\mathrm{2}}^{z_H+\mathrm{\Delta }z/\mathrm{2}};\end{array}$$

this follows the numerical solution to within ∼20 % for $\mathrm{0.3}\lesssim z/h_{\mathrm{mE}}\lesssim \mathrm{2}$, more so 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−φ₀). 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 Eqs. (25) and (27), while the latter is obtained via Eq. (29).

One can see in the left-hand plot of Fig. 1 that the Ekman solution produces effectively less mixing and consequently a higher shear exponent than Ellison's. Similarly, away from the surface (for $z>{\mathit{\nu }}_{\mathrm{Ek}}/\left(\mathit{\kappa }{u}_{*}\right)$, i.e., $z/h_{\mathrm{Ek}}\gtrsim \mathrm{0.1}$) in the right-hand plots of Fig. 1 one can see the dimensionless Ekman veer exceeds that predicted by the Ellison solution; this is consistent with the Ekman ABL turning angle of ${\mathit{\gamma }}_{\mathrm{0},\mathrm{Ek}}=\mathrm{45}$^∘, which exceeds the γ₀ of 5–15^∘ predicted by Ellison's form (van der Laan et al., 2020). However, we note that the depth h can differ between the Ekman and Ellison solutions; h_Ek=h_mE only if one chooses ${\mathit{\nu }}_{\mathrm{Ek}}=(\mathit{\kappa }{u}_{*}{)}^{\mathrm{2}}/\mathrm{2}f$. We also point out that for larger c_G (not shown in figure), near the surface ($z/h_{\mathrm{mE}}\lesssim \mathrm{0.1}$) Ellison's veer grows 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.

2.4 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 crosswind stress 〈vw〉_⟂ appear in Eqs. (14) and (16), along with $\left|G\right|$. In wind energy applications, engineers typically lack site-specific stress profiles unless they are taken from flow modeling; 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 $\left|S\right|/\left|G\right|$ in Eqs. (14)–(16). To obtain a practical form relating shear to veer, we can start by parameterizing $\left|S\right|/\left|G\right|$; fortunately $\left|G\right|$ 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

$$\begin{array}{}\text{(31)}& \left|G\right|={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\sqrt{{\left[\ln \left({\displaystyle \frac{u_{*}}{f{z}_{\mathrm{0}}}}\right)-A\right]}^{\mathrm{2}}+B^{\mathrm{2}}},\end{array}$$

with components

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\sin \left({\mathit{\phi }}_{\mathrm{G}}-{\mathit{\phi }}_{\mathrm{0}}\right)& =-B{\displaystyle \frac{u_{*}}{\mathit{\kappa }\left|G\right|}}\mathrm{and}\cos \left({\mathit{\phi }}_{\mathrm{G}}-{\mathit{\phi }}_{\mathrm{0}}\right)\\ & ={\displaystyle \frac{u_{*}}{\mathit{\kappa }\left|G\right|}}\left[\ln \left({\displaystyle \frac{u_{*}/f}{z_{\mathrm{0}}}}\right)-A\right],\end{array}\end{array}$$

where the empirical coefficients {A,B} are assumed to be constants in typical wind application. The geostrophic drag coefficient $c_{\mathrm{G}}\equiv u_{*}/\left|G\right|$ and ABL turning (cross-isobar angle) φ_G are seen to vary with surface Rossby number Ro₀ (Blackadar and Tennekes, 1968); these and {A,B} have been shown to depend on dimensionless stability $L^{-\mathrm{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 the influence of the capping inversion strength.⁹ With such assumptions, one can also write an (approximate) “reverse” form of Eq. (31) to get the drag coefficient as (Troen and Petersen, 1989)

$$\begin{array}{}\text{(33)}& c_{\mathrm{G}}\simeq {\displaystyle \frac{c_{\mathrm{rGDL}}}{\ln R{o}_{\mathrm{0}}-A}},\end{array}$$

where the surface Rossby number is $R{o}_{\mathrm{0}}\equiv \left|G\right|/\left(f{z}_{\mathrm{0}}\right)$ and c_rGDL is taken to be 0.485 following its use in the wind resource program WAsP for several decades. Alternate forms of Eq. (33) exist, such as that of Hess and Garratt (2002); the latter corresponds simply to setting A=1.28 and c_rGDL=0.472 in Eq. (33). For a given roughness length z₀ and measured wind speed $\left|S\right|$, lacking the (surface) friction velocity u_*, one needs a relation to connect u_* and $\left|S\right|$ in order to get $\left|G\right|$. 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_{*}=\mathit{\kappa }\left|S\right|/\ln (z/z_{\mathrm{0}})$ within Eq. (31) or alternately $\left|S\right|/\left|G\right|=(c_{\mathrm{G}}/\mathit{\kappa })\ln (z/z_{\mathrm{0}})$ using Eq. (33), where in the latter Eq. (31) is also employed to find $\left|G\right|$ within Ro₀.

In practice one would like a direct estimate for the veer, using the routinely measured shear, since α is seen to drive $\partial \mathit{\phi }/\partial z$. One way could be to just ignore the stress divergence terms in Eqs. (14) or (16), which with calculation of $\left|G\right|$ mentioned just above considerably simplifies the problem. However, this might not be justified, particularly if $u_{*}^{\mathrm{2}}/\left(fh\right)$ is not negligible compared to $\left|S\right|$, as seen from comparing contributions to Eqs. (14)–(16); this can be seen using the scaling $\partial \langle uw\rangle /\partial z\approx u_{*}^{\mathrm{2}}/h$, where h is the ABL depth (e.g., Wyngaard, 2010). Thus we consider estimating vertical derivatives of the stresses, starting with the $\partial \langle uw\rangle /\partial z$ just mentioned, which can be used in Eq. (16). Similarly, one can estimate ${\partial }^{\mathrm{2}}\langle vw\rangle /\partial z^{\mathrm{2}}\approx c_{vw}{u}_{*}^{\mathrm{2}}/h^{\mathrm{2}}$ or

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}\langle vw\rangle /\partial z^{\mathrm{2}}}{f\left|G\right|}}\approx c_{vw}{\displaystyle \frac{u_{*}^{\mathrm{2}}}{f\left|G\right|h^{\mathrm{2}}}}=c_{vw}{\displaystyle \frac{c_G^{\mathrm{2}}}{h}}R{o}_h,\end{array}$$

where $R{o}_h\equiv G/\left(fh\right)$ is the Rossby number based on ABL depth and of order 1; we will treat c_vw as an empirical constant which is tuned later below. To use Eq. (16) we also need to find sin φ_G; employing Eq. (32) and using trigonometric identities to expand sin (φ−φ₀), with some rearrangement one obtains

$$\begin{array}{}\text{(35)}& \sin {\mathit{\phi }}_{\mathrm{G}}={\displaystyle \frac{c_{\mathrm{G}}}{\mathit{\kappa }}}\left\{\left[\ln \left({\displaystyle \frac{u_{*}}{f{z}_{\mathrm{0}}}}\right)-A\right]\sin {\mathit{\phi }}_{\mathrm{0}}-B\cos {\mathit{\phi }}_{\mathrm{0}}\right\}.\end{array}$$

Employing this, Eq. (34), and $\partial \langle uw\rangle /\partial z\approx u_{*}^{\mathrm{2}}/h$, along with Eqs. (31) or (33), allows one to then use Eq. (16).

On the other hand, using Eq. (14) is simpler and more convenient than Eq. (16) because it only requires $\partial \langle vw\rangle /\partial z$ in addition to the second derivative of 〈vw〉 just approximated in Eq. (34) above, so one can also simply approximate $\partial \langle vw\rangle /\partial z\approx h{\partial }^{\mathrm{2}}\langle vw\rangle /\partial z^{\mathrm{2}}$ and use Eq. (34); the GDL forms Eqs. (31) and (33) which then allow one to get Ro_h and c_G, respectively. Whether using Eq. (16) or (14), we note that the shear contribution to veer includes a surface Rossby number (Ro₀) dependence through $S/\left|G\right|$, 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.

3 Analysis and discussion

This section presents an analysis of results from RANS simulations of the neutral atmospheric boundary layer¹⁰ 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₀ and Ro_h) encountered in nature. An investigation of observations, spanning turbine rotor heights for six different flow regimes and conditions across four locations, 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.

Figure 2 Demonstration of one-dimensional RANS solver results, conforming to Eq. (13) (a) and Eq. (9) (b). Dashed line represents 1:1 prediction; simulated ABL depth h_eff calculated from Eq. (36).

3.1 RANS simulations of neutral ABLs

3.1.1 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 computational fluid dynamics (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 modeled via the length-scale limiter model of Apsley and Castro (1997) implemented into the k–ε turbulence closure equations solved by Ellipsys1D, as outlined in van der Laan et al. (2020); this includes the 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 occurring in the RANS equations, via the flux-gradient relation and ${\mathit{\nu }}_{\mathrm{T}}=C_{\mathit{\mu }}{k}^{\mathrm{2}}/\mathit{\epsilon }$; thus we see that such turbulence closure gives stresses aligned with velocity gradients.

The domain height is set to 10⁵ 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 ($\mathrm{d}k/\mathrm{d}z=\mathrm{0}$) and ε is set to the logarithmic value, and 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 $\left|G\right|$ at all heights; the ABL depth grows upward until convergence occurs, providing a steady solution and h for a given choice of z₀, pressure gradient (thus G and f), and turbulence (k–ε) limiting length scale ℓ_max. The Buckingham pi theorem can be used to reduce the four parameters $\mathit{\{}{z}_{\mathrm{0}},G,f,{\mathrm{\ell }}_{\max }\mathit{\}}$ into two dimensionless groups, namely Rossby numbers for z₀ and ℓ_max; for length-scale-limited k–ε RANS in the neutral ABL, one further has the relation (van der Laan et al., 2020)

$$\begin{array}{}\text{(36)}& h={\mathrm{\ell }}_{\max }^{\mathrm{0.6}}{\left({\displaystyle \frac{\left|G\right|}{\left|f\right|}}\right)}^{\mathrm{0.4}},\end{array}$$

thus giving us the two Rossby numbers Ro₀ 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${}^{\mathrm{5}}<R{o}_{\mathrm{0}}<{\mathrm{10}}^{\mathrm{10}}$ and $\mathrm{15.8}<R{o}_h<\mathrm{661}$. For simplicity $\left|G\right|$ was set to 10 m s⁻¹ and f to 10⁻⁴ s⁻¹ in the simulation set spanning these ranges of Rossby numbers. However, note that Rossby similarity means that for a given pair of {Ro₀,Ro_h} and {z₀,h} one has many (infinite) combinations of $\mathit{\left\{}\right|G|,f\mathit{\}}$ which give the same $\left|G\right|/f$ and thus the same dimensionless profile shapes of velocity, i.e., speed and direction as a function of dimensionless height $zf/\left|G\right|$. At any rate, the simulations cover ranges of (exceeding) the following: ABL depths of 200–2000 m, roughness lengths from water's roughness (0.1 mm) up to 2.5 m, and $\left|G\right|$ from 5–50 m s⁻¹.

3.1.2 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 Eq. (7). Figure 2 displays both sides of Eqs. (9) and (13) for four cases representing somewhat common real-world conditions for heights between 50–200 m.

From Fig. 2 one can see that the Ellipsys1D solutions conform to Eqs. (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. Admitting 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 to 150 m. Figure 3 displays $\mathrm{\Delta }\mathit{\phi }{|}_{\mathrm{50}\mathrm{m}}^{\mathrm{150}\mathrm{m}}$ plotted over different roughnesses for the two ABL depths represented in the cases shown in the previous figure, namely 490 and 1047 m.

Figure 3 Roughness dependence of turning (in degrees clockwise) seen for two representative ABL depths from one-dimensional RANS simulations over a range of roughness lengths plotted directly against z₀ (a) and alternately versus $\mathrm{1}/\ln \left(R{o}_{\mathrm{0}}\right)$ (b). Lines in (b) indicate linear trend.

From the right-hand plot in Fig. 3 one can see that Δφ is roughly proportional to $\mathrm{1}/\ln \left(R{o}_{\mathrm{0}}\right)$, as expected from the $S/\left|G\right|$ contribution to veer considering Eqs. (31) and (33).

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

Looking back on Eq. (34), we may also expect a Ro_h dependence in the veer, at least considering the stress gradient contributions. Figure 4 shows veer across three different rotor extents (z=50–100, 50–150, 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 only for one roughness because the curves of veer versus ABL depth and Ro_h look nearly identical when using any other z₀ (or Ro₀) value, such as water roughnesses less than 0.3 mm. In other words, the sensitivity of veer to h is essentially independent of z₀ if one varies these separately from case to case as in our numerical simulations. Looking at these results, we note a behavior that is consistent with the estimates for stress-gradient contributions following Eq. (34): as indicated by the dotted lines in Fig. 4, the veer is empirically found to be proportional to $R{o}_h^{\mathrm{1.4}}$ (or h^−1.4) over a range of ABL depths routinely observed in reality (h∼200–800 m); the dependence softens to be linear in Ro_h (or $\mathrm{1}/h$) for depths approaching h∼1 km, which are also commonly observed in nature (e.g., Liu and Liang, 2010). For yet deeper ABLs which are more rarely encountered, the height dependence vanishes; this can be intuitively interpreted, as $\mathrm{\Delta }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₀: across all Ro₀ a variation in $\mathrm{\Delta }\mathit{\phi }{|}_{\mathrm{50}\mathrm{m}}^{\mathrm{150}\mathrm{m}}$ of only several degrees is seen, whereas across the common range of Ro_h a variation of more than 15^∘ is shown.

Figure 5 Profiles of contributions to $\mathrm{d}\cos \mathit{\gamma }/\mathrm{d}z$ in Eq. (12) due to shear (a), stress gradients with Coriolis (b), and their sum (c). Five RANS simulations shown (two roughnesses and two ABL depths over land, one over sea) over typical turbine rotor heights; the listed z₀ and h correspond to Rossby numbers using G=10 m s⁻¹ and $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹.

Now that we have seen in Figs. 3 and 4 how the veer (or simply the turning Δφ for typical rotor Δz) depends on z₀ and h, presumably due to the $S/\left|G\right|$ (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 in cos γ, as indicated by Eqs. (12) and (13). Accordingly, Fig. 5 presents the two contributions to the dimensionless veer $\partial \cos \mathit{\gamma }/\partial z$ derived in Eq. (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.

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. (a) Using G=10 m s⁻¹ and $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ over range of ABL depths spanning the values used in Figs. 2–5 for water and typical land roughness. (b) Again with G=10 m s⁻¹ and $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ over range of z₀ spanning those used in Figs. 2–5 for two ABL (typical) depths used in previous figures. (c) Over wider range of {Ro₀,Ro_h} spanning that found in nature; note larger vertical axis scale.

One can note from Fig. 5 that the shear and stress-gradient and Coriolis contributions largely offset each other, with each being an order of magnitude larger than their sum, which is equal to the dimensionless veer $\partial \cos \mathit{\gamma }/\partial z$. The vertical profiles of “point-wise” veer shown in the figure, which were calculated using third-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₀(z₀) having a smaller impact than Ro_h(h).

This can be put into a more practical context by considering the variation in shear and veer together across the range of Rossby numbers found in atmospheric flows. Figure 6 displays turning versus shear exponent, with each calculated across Δz from 50–150 m. The figure shows three plots of {α,Δφ}: one for a range of Ro_h equivalent to h ranging from 490 to 1047 m over two different z₀ (land and sea), one for a range of Ro₀ equivalent to z₀ values varying from 0.016–25 cm for two different ABL depths h (which bracket the range of h in the left-hand plot), and one over the entire atmospheric range of both Ro₀ and Ro_h.

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₀ 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 subsection, the right-hand plot is more in line with observations despite the RANS solutions representing nominally neutral conditions¹¹ over uniform surfaces with neglect of shear-stress misalignment and baroclinity.

Figure 7 (a) PDF of shear exponent α, weighted by frequency of occurrence for different atmospheric stability conditions; neutral is defined by $\left|L^{-\mathrm{1}}\right|<\mathrm{0.001}$ m⁻¹, stable has $L^{-\mathrm{1}}>\mathrm{0.001}$ m⁻¹, and unstable has $L^{-\mathrm{1}}<-\mathrm{0.001}$ m⁻¹. (b) Joint PDF of shear exponent calculated from 60–160 m height and inverse Obukhov length (surface-layer stability) at z=10 m from sonic anemometers over the homogeneous land sectors at Høvsøre, along with M–O similarity implied by Eq. (2) shown by dashed cyan line; JPDF value of 1000 corresponds to approximately 3.5 % occurrence. Measurements span one decade, starting in 2005.

3.2 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 effects of stability. The datasets are the same as those 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 (Peña, 2019) for two virtual rotor spans, from 45–140 and 80–200 m over 1 year; the Dutch research site Cabauw (Beljaars and Bosveld, 1997), from 80–200 m height for 2 years; and 1 year from a commercial site dubbed “MR” which sits on a ridge over a mostly forested ($>\sim \mathrm{3}/\mathrm{4}$) area but dominated by hills having elevation differences up to ∼200 m within 10 km distance, using anemometers at 40–136 m height.¹²

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.

3.2.1 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⁻¹) 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 Eqs. (4) and (5), α tends to correlate with stability ($\mathrm{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 $\mathit{\alpha }{|}_{\mathrm{60}\mathrm{m}}^{\mathrm{160}\mathrm{m}}$ and L⁻¹ 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 min averages over a 10-year period.

From Fig. 7 one sees the cloud of observed $\mathit{\{}\mathit{\alpha },L^{-\mathrm{1}}\mathit{\}}$ follow somewhat the curve of ${\mathrm{\Phi }}_m/\left[\ln \right(z/z_{\mathrm{0}})-{\mathit{\psi }}_m(z/L\left)\right]$ implied by M–O theory and Eq. (2) but with most α exceeding the similarity-based form; the shear exceeds M–O theory's prediction primarily due to the upper-level height (160 m) being above the surface layer.¹³ The left-hand panel of Fig. 7 also shows the distribution of α for neutral ($\left|L^{-\mathrm{1}}\right|<\mathrm{0.001}$ m⁻¹), stable ($L^{-\mathrm{1}}>\mathrm{0.001}$ m⁻¹), and unstable ($L^{-\mathrm{1}}<-\mathrm{0.001}$ m⁻¹) flow regimes, weighted by frequency of occurrence to show the relative contributions to the overall distribution. The threshold of ±0.001 m⁻¹ for L⁻¹ is a sensible choice because then $z/\left|L\right|\ll \mathrm{1}$ (consistent with neutral conditions) in the surface layer, which is generally taken to have a thickness of 100 m or less (roughly $h/\mathrm{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 more so in unstable and yet less often in neutral conditions; overall, α<0 occurs less than 5 % of the time over the 60–160 m span 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).

3.2.2 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 homogeneous offshore/open-fetch ($\mathrm{240}{}^{\circ }<\mathit{\phi }<\mathrm{300}{}^{\circ }$) and over-land ($\mathrm{60}{}^{\circ }\mathit{\phi }<\mathrm{120}{}^{\circ }$) directions. Shear and veer are shown calculated over height spans of 60–160 m as well as 100–160 m in the figure, which is provided to show the statistical and behavioral differences between shear and veer.

Figure 8 Distributions of shear exponent (a) and corresponding veer (b) 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 and dashed for those spanning 100–160 m.

From the two plots in Fig. 8 one can see that the most common α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$, i.e., the portions of P(α) and $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ between these two height ranges. This happens over both land and sea, though both α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ vary with height more for the offshore flow than for the homogeneous land directions. The change in mean shear exponent 〈α〉 from 60–160 to 100–160 m is less than +5 % over land, and $<+\mathrm{30}$ % is seen over sea, while the mean veer $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\rangle$ is seen to increase by factors of $\sim \mathrm{5}/\mathrm{3}$ and 2 over land and sea, respectively.

Figure 9 (a) PDF of veer $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$, weighted by frequency of occurrence for different stability conditions; neutral has $\left|L^{-\mathrm{1}}\right|<\mathrm{0.001}$ m⁻¹, stable has $L^{-\mathrm{1}}>\mathrm{0.001}$ m⁻¹, and unstable has $L^{-\mathrm{1}}<-\mathrm{0.001}$ m⁻¹. (b) Joint probability distribution of veer calculated from 60–160 m height and inverse Obukhov length (surface-layer stability) at z=10 m from sonic anemometers over the homogeneous land sectors at Høvsøre; here a JPDF value of 3000 corresponds to 4.2 % occurrence. Measurements span one decade, starting in 2005.

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(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ peak does not differ significantly from land to sea here. Further, the (logarithmic) slope of $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ versus $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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₀ primarily affects α, while $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 cause 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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ are larger for the 100–160 m measurements than those from 60–160 m; this is again due to the greater 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(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ for neutral ($\left|L^{-\mathrm{1}}\right|<\mathrm{0.001}$ m⁻¹), stable ($L^{-\mathrm{1}}>\mathrm{0.001}$ m⁻¹), and unstable ($L^{-\mathrm{1}}<-\mathrm{0.001}$ m⁻¹) flow regimes, weighted by frequency of occurrence (indicating relative contributions to the full PDF), as well as the joint distribution of stability and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$.

From Fig. 9 one sees that in comparison with P(α) shown in Fig. 7, the peaks of veer distributions $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ do not depend so much on stability. However, as with the shear distribution, $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 is 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(\mathit{\alpha },\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ one sees that for the most common veer values ($\mathrm{0}\lesssim \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\lesssim \mathrm{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.

To show the behavior of veer across different locations, Fig. 10 displays the PDFs of veer from a number of sites, all of which have similar Δz and cover typical turbine rotor extents.

Figure 10 Probability density function (distribution) of veer, $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$, for all conditions at the various sites and cases considered.

From Fig. 10 we see that for veer magnitudes exceeding the 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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\gtrsim \mathrm{0.2}$^∘ m⁻¹ in each case. These slopes correspond to (conditional) PDFs for the largest veer of the form

$$\begin{array}{}\text{(37)}& P\left({\displaystyle \frac{\mathrm{\Delta }\mathit{\phi }}{\mathrm{\Delta }z}}|{\displaystyle \frac{\mathrm{\Delta }\mathit{\phi }}{\mathrm{\Delta }z}}>\mathrm{mode}\left\{{\displaystyle \frac{\mathrm{\Delta }\mathit{\phi }}{\mathrm{\Delta }z}}\right\}\right)\propto \exp \left[{\displaystyle \frac{-\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z}{{\mathrm{{\rm Y}}}_{\mathrm{veer}}}}\right],\end{array}$$

where the characteristic veer scale defined by ${\mathrm{{\rm Y}}}_{\mathrm{veer}}^{-\mathrm{1}}\equiv \partial \left[\ln P\right(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\left)\right]/\partial (\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ ranges from roughly 0.07 to 0.11^∘ m⁻¹. The lowest Υ_veer corresponds to the offshore Høvsøre case, while the highest Υ_veer matches the Østerild case from 45–140 m. We expect larger Υ_veer to correspond to occurrences of higher $\mathrm{1}/L$, i.e., a larger width σ₊ of the stable-side distribution $P(\mathrm{1}/L)$ following Kelly and Gryning (2010); essentially the large-veer PDF in Eq. (37) is conditional on stable conditions; i.e., we could express it as $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|L^{-\mathrm{1}}>\mathrm{0})\propto \exp [-(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)/{\mathrm{{\rm Y}}}_{\mathrm{veer}}]$. The dominance of stable conditions reported by Peña (2019) for z≳100 m at Østerild is consistent with this, though the data from z=80–200 m (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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\gtrsim \mathrm{0.4}$^∘ m⁻¹, there is consistency: the two largest Υ_veer values occur for z=80–200 and z=45–140 m. 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(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ in Fig. 10 are at smaller $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ for the forest-dominated Østerild cases, with the peak of the offshore Høvsøre veer distribution falling between these and the $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ corresponding to the land cases of Høvsøre, Cabauw, and MR. Note that the most commonly found veer values are generally dominated by neutral conditions (or modestly stable for the exceptional Østerild site above 100 m) and point out that the mode of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ is essentially the same (0.005–0.006^∘ m⁻¹) for the land cases that are not dominated by forest. Further considering the RANS simulation results from Fig. 6 discussed earlier, the mode of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 modes of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 2–5 km 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 min mean wind speeds, $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z,U)$. Along with the joint distribution, the mean veer conditioned on wind speed, $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\rangle |U$, is displayed.

Figure 11 Joint distribution of veer and wind speed at sites considered. Solid line shows $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\rangle |U$, calculated using 1 m s⁻¹ bins; lightest shades are 2 % as likely as darkest color in each plot.

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⁻¹ for the offshore case) because sea–air heat fluxes and associated $\mathrm{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 pre-construction annual energy production (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 in 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 and cases considered. It also adds a line to show the over-land Høvsøre case filtered for neutral conditions.

Figure 12 Measured standard deviation of veer conditioned on wind speed, again using 1 m s⁻¹ bins, for the sites considered. Neutral conditions at Høvsøre are defined as in earlier figures, i.e., $\left|L^{-\mathrm{1}}\right|<\mathrm{0.001}$ m⁻¹.

Consistent with the joint-PDFs $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 more so over land than water. It is also seen that for the onshore Høvsøre case ${\mathit{\sigma }}_{(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z}$ and mean wind speed 〈U〉, the rms veer conditioned on wind speed roughly follows the empirical form

$$\begin{array}{}\text{(38)}& {\mathit{\sigma }}_{(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)|U}={\left[〈{\left({\displaystyle \frac{\mathrm{\Delta }\mathit{\phi }}{\mathrm{\Delta }z}}\right)}^{\mathrm{2}}|U〉\right]}^{\mathrm{1}/\mathrm{2}}\approx {\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z}\exp \left[{\displaystyle \frac{-U}{\langle U\rangle }}\right]\end{array}$$

up to about 12 m s⁻¹ over land and to higher speeds offshore. A more complicated speed-dependent variability in veer is seen for the MR case, with higher ${\mathit{\sigma }}_{(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)|U}$ at speeds above 15 m s⁻¹ caused (presumably) by hill-induced turning. This has two consequences worth mentioning: firstly, that turbines at a site such as MR can experience persistent veer above rated speed, potentially increasing loads and/or reducing power below rated; and 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.

3.3 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 and 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 Sect. 2.3 and 2.4.

Figure 13 Top panels: 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 occurrence rate normalized per wind speed range (each plot has a different color scale, showing occurrence rate in increments of $\mathrm{1}/\mathrm{10}$, with the lightest shade representing 10 % as likely as the darkest shade). Bottom panels: the same joint distributions shown with unscaled rate of occurrence (number of counts per $\mathit{\{}\mathit{\alpha },\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\mathit{\}}$ bin, 2000 corresponds to about 3.9 %); axis ranges are chosen to compare with later figures.

Following the previous subsection, we first consider the joint behavior of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z,\mathit{\alpha })$ in neutral conditions, over typical turbine operation speeds (4–25 m s⁻¹), and separately over different speed ranges (4–8, 8–12, 12–16, and 16–25 m s⁻¹); counts are used instead of PDF per wind speed range to show relative frequencies of occurrence.

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 the bottom of Fig. 13. All plots use the same color scale; color bar denotes count, where 3000 corresponds to about 1.4 % occurrence.

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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ and α for U>16 m s⁻¹; 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 nonstationarity and particularly sampling uncertainty; per the latter the integral timescale increases roughly as U⁻¹ (Wyngaard, 2010), so fewer integral timescales are “sampled” per each 10 min period. This is also evident considering the previous plot of ${\mathit{\sigma }}_{(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)|U}$ versus U in Fig. 12, where one sees ${\mathit{\sigma }}_{(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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⁻¹. Also, the overall JPDF (joint PDF) $P(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z,\mathit{\alpha })$ appears similar to that in the most common speed range (8–12 m s⁻¹). Aside from nonstationarity and sampling effects one does not expect much speed dependence in neutral conditions, considering the α-related part of Eqs. (14)–(16) behaves as $\left|S\right|/G$, which following Eq. (33) has a weak $\left|S\right|$ dependence through $(\ln R{o}_{\mathrm{0}}-A{)}^{-\mathrm{1}}$; the RANS results also confirm this. We note a joint trend between α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 15 Joint distribution P(Δφ,α) at sites considered. Solid lines: mean veer conditioned on shear exponent, $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$; dotted lines: simple estimate via shear portion using Eq. (39); dashed lines: estimate including estimate of crosswind stress–Coriolis contribution, Eq. (40). Lightest shades are 10 % as likely as darkest shades in each plot.

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 occurrence of higher veer and shear, as well as negative α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ with α for larger α due to stable conditions. One can see that at the most common (8–12 m s⁻¹) and lower wind speeds, which occur in the range below rated speed for typical turbines, there is a significant increase in $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ with α in the more stable conditions where α≳0.3; this higher “slope” of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ versus α is likely enhanced by the shallower ABLs which generally occur along with stable surface-layer conditions (note that the stability metric L⁻¹ was measured in the ASL), whereby additional 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 $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z=\mathrm{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., $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$, as solid lines).

From this figure we see a number of trends across the six cases analyzed. First, some nonlinear variation in veer with α is evident, along with the (less common) occurrence 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., $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ exceeds the most commonly observed values of $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$; however, the MR site does not show such skewed behavior (consistent with Fig. 10), presumably due to the complex terrain there. We note the conditional mean veer $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ is also more clearly nonlinear in α, becoming less dependent on α in low-shear (and negative) conditions; the MR site is an exception to this, with hill-induced height-dependent turning causing larger veer for α smaller than the most commonly observed values there.

3.3.1 Simplified estimate of veer per α

Figure 15 also includes two predictions based on the theory presented earlier. First, as discussed at the end of Sect. 2.4, using only the shear-associated ($\left|S\right|/\left|G\right|$) portion of Eq. (14) to be practical, we arrive at the estimate

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}\approx {\displaystyle \frac{\left|S\right|}{\left|G\right|}}{\displaystyle \frac{\mathit{\alpha }}{z}}/\sqrt{\mathrm{1}-{\left[{\displaystyle \frac{\left|S\right|}{\left|G\right|}}\right]}^{\mathrm{2}}},{\displaystyle \frac{\left|S\right|}{\left|G\right|}}\approx c_{s\mathit{\alpha }}{\displaystyle \frac{c_{\mathrm{rG}}}{\mathit{\kappa }}}{\displaystyle \frac{\ln \left(z/z_{\mathrm{0}}\right)}{\left(\ln R{o}_{\mathrm{0}}-A\right)}};\end{array}$$

compared to Eq. (14) the negative sign has been dropped to express the veer in coordinates commonly used in wind energy, i.e., clockwise positive. The basis for the simple shear-driven form can be understood by recalling Sect. 3.1 and Fig. 5, where we showed that the shear and crosswind stress-curvature contributions behaved in nearly identical but opposite fashions, with their sum amounting to $\mathrm{d}\cos \mathit{\gamma }/\mathrm{d}z$; Eq. (39) can be considered a simple model assuming the veer behaves like either of its two components but is simply smaller in magnitude. The practical form of $\left|S\right|/\left|G\right|$ in Eq. (39) employs the log-law for wind profile and reverse geostrophic drag law Eq. (33) for $u_{*}/\left|G\right|$. The constant c_sα crudely accounts for the (competing) effects of stability on both $\left|S\right|$ and the geostrophic drag (and any other mechanisms affecting $\left|S\right|/\left|G\right|$) but also accounts for the smaller magnitude of $\partial \mathit{\phi }/\partial z$ compared to its shear-driven component. Within the surface Rossby number Ro₀, the geostrophic speed G is calculated using Eq. (32), wherein u_* is found via the log-law and $\left|S\right|$ with z₀. To make the plots of Eq. (39) in Fig. 15 for each site, the $\left|S\right|$ is calculated per each bin of α, with the case-specific parameters $\mathit{\{}{z}_{\mathrm{0}},f,z\mathit{\}}$ used as well. At any rate, the practical parameterization using c_sα with the log-law and (neutral) reverse GDL in Eq. (39) can roughly fit the mean conditional veer at and 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₀=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 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 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 (labeled σ+ following Kelly and Gryning, 2010); we 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⁻¹ values consistent with observed distributions) 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 conditions (not shown), in contrast to the value of 0.7 which fits when all stabilities are considered there.

Figure 16 Statistics of veer conditioned on shear exponent across sites considered.

3.3.2 Veer estimate including both α and crosswind stress

Note that for simplicity, Eq. (39) ignored the effect of crosswind stress; it neglects not only 〈vw〉 but consequently also Ro_h, though it does incorporate the effect of Ro₀ seen in the simulations of Sect. 3.1. Thus we also consider an approximation of the 〈vw〉 terms using Eq. (34) in Eq. (14), which introduces Ro_h, along with the parameterization for $\left|S\right|/\left|G\right|$ from Eq. (39):

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial z}}\approx {\displaystyle \frac{\frac{\left|S\right|}{\left|G\right|}\frac{\mathit{\alpha }}{z}+c_{vw}\frac{c_{\mathrm{G}}^{\mathrm{2}}}{h}R{o}_h}{\sqrt{\mathrm{1}-{\left[\frac{\left|S\right|}{\left|G\right|}+c_{vw}{c}_{\mathrm{G}}^{\mathrm{2}}R{o}_h\right]}^{\mathrm{2}}}-c_{\mathrm{G}}^{\mathrm{2}}R{o}_h}}\\ & {\displaystyle \frac{\left|S\right|}{\left|G\right|}}\approx c_{s\mathit{\alpha }}^{\prime }{\displaystyle \frac{c_{\mathrm{rG}}}{\mathit{\kappa }}}{\displaystyle \frac{\ln \left(z/z_{\mathrm{0}}\right)}{\left(\ln R{o}_{\mathrm{0}}-A\right)}},\end{array}\end{array}$$

where c_G is found using Eq. (33), $\left|S\right|/\left|G\right|$ is calculated the same way as done earlier for Eq. (39), and $\left|G\right|$ within Ro_h is calculated as it was within Ro₀ of Eq. (39). To use Eq. (40) the ABL depth must be prescribed, along with the constant c_vw and the parameters $\mathit{\{}z,|S|,z_{\mathrm{0}},f\mathit{\}}$ also employed for Eq. (39). Given the negative curvature of lateral stress, ${\partial }^{\mathrm{2}}\langle vw\rangle /\partial z^{\mathrm{2}}<\mathrm{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 Eq. (39). 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\mathit{\alpha }}^{\prime }$ within $\left|S\right|/\left|G\right|$ is slightly larger than c_sα of Eq. (39) in order for Eq. (40) to fit the observed $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$; the values of c_sα=0.5 are replaced by $c_{s\mathit{\alpha }}^{\prime }=\mathrm{0.7}$, and c_sα=0.7 and 0.8 for Høvsøre and Cabauw are replaced by $c_{s\mathit{\alpha }}^{\prime }$ of 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 and 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 $\langle \partial \mathit{\phi }/\partial z|\mathit{\alpha }\rangle$ using Eq. (40) are not better than the simpler form of Eq. (39), 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 wish to simply predict veer based on α; to this end, for practical applicability we suggest using Eq. (39). 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 in Ro_h or Ro₀ at a given site. The spread (scatter) around the mean veer seen in Fig. 15 is due to variation in stability as well as Ro_h or Ro₀, and variation from site to site is also due to different distributions of Ro_h or Ro₀; 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 and cases considered.

One immediately sees the character of $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ 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 note 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 ($\mathit{\alpha }<\sim \mathrm{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 that the long-term variability in veer ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }}$ is lower offshore for the most commonly occurring α there, while veer variability does not differ so much for the most common conditions across the other sites and 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 ($\mathit{\alpha }>\sim \mathrm{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 ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }}/\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ 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 occurring shear; the more homogeneous sites and cases exhibit comparable ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }}/\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ 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 and 16 actually incorporated mean speed conditioned on α, though use of each site's corresponding overall mean speed $\langle \left|S\right|\rangle$ gave nearly identical results as those shown in the plots (within 2 %, not shown).

4 Summary and conclusions

We have derived relationships between shear exponent (α) and veer ($\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$) in a manner which avoids atmospheric-surface-layer (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 roughness z₀ and ABL depth h is also elucidated (through Rossby numbers Ro₀ and Ro_h defined by each) through numerical solution of the one-dimensional RANS equations under neutral conditions with length-scale-limited k–ε turbulence closure (i.e., neutral but also translatable to stable conditions; see van der Laan et al., 2020).

The derived equations and RANS results essentially show that veer most simply arises from two contributions: the shear and the vertical variation in crosswind shear stress at a given height (mostly through ${\partial }^{\mathrm{2}}\langle vw\rangle /\partial z^{\mathrm{2}}$ but also via $\partial \langle vw\rangle /\partial 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 $\mathrm{\Delta }z/h$; in contrast, $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ more strongly depends on the ABL depth h, increasing as $R{o}_h^n$, where n is between 1 and 1.4 for the h most commonly encountered in nature (though $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ does also vary with $\mathrm{1}/\ln R{o}_{\mathrm{0}}$). These behaviors are consistent with the shear–veer relations derived in Sect. 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; note that RANS solutions using mixing-length-type closures (as well as, e.g., Weather Research and Forecasting planetary boundary layer (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 crosswind stress.

The actual “real-world” behavior of shear exponent and veer has also been investigated from multi-year measurements at four sites covering 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 $\mathit{\{}\mathit{\alpha },\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z\mathit{\}}$ 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(\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$, while stable conditions are responsible for large α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$; neutral conditions contributed mostly to the peaks of the shear and veer distributions. Stability effects are consequently seen to increase the long-term variability in veer and shear, as well as veer for a given α – particularly for the commonly occurring wind speeds which tend to occur below the rated speed of modern wind turbines (e.g., Appendix B of Kelly and Jørgensen, 2017). The mean of both α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ was larger compared to neutral conditions due to stably stratified conditions enhancing α and $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z$ more than unstable conditions (we note that sites having a distribution of $\mathrm{1}/L$ more dominated by unstable conditions, possibly some offshore, could have mean behavior similar to that found in neutral conditions).

A 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 with the capping inversion, $\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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).

Two practical shear–veer relationships were derived, including parameterizations for typically unmeasured quantities contained within them, and then compared to the joint distributions $P(\mathit{\alpha },\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z)$ and the $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ measured from all sites over all conditions. A simplified form Eq. (39) neglecting the stress contributions was tested, as well as one Eq. (40) containing the crosswind 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 $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$: basically one coefficient in Eq. (39) and two in Eq. (40), all of which were of order 1 and universal (constant) across all six sites and flow situations analyzed. The form Eq. (40) for veer including crosswind stress did not give a better match to observations of $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ across sites, compared to the simpler formula (Eq. 39), 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) for $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ during the most commonly observed (moderate speeds and shear) and highest-impact (large-veer and stable) situations, with underpredictions of mean veer occurring in low-veer conditions. The observed $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ is nonlinear in α, whereas the derived forms were nearly linear, with the inclusion of crosswind 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 underpredictions for smaller α, dominated by unstable conditions, evoke this. Consistent with this, the hilly MR site shows yet more low-shear deviation from our predictions due to inhomogeneity-related horizontal transport (recalling that low shear means less shear production of TKE).

Beyond the comparison of derived analytical forms with measurements of conditional mean veer $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$, some general trends were also noted. For a given α≳0.2, $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ was larger offshore than for the onshore cases (though note that larger α is 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 $\langle \mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }\rangle$ was smaller than other sites. As for the mean veer, for α≳0.3 the long-term variability ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }z|\mathit{\alpha }}$ was also found to be larger offshore; this may have an 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 (${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\phi }/\mathrm{\Delta }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 stability within $\left|S\right|/\left|G\right|$ 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 and within both the shear and cross-stress contributions in concert with the stability-perturbed geostrophic drag law (Arya and Wyngaard, 1975; Kelly and Troen, 2016). Future work includes the incorporation of 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 occurring 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 load calculations.