Atmospheric stability is our friend (yes it is)

When dealing with site condition studies for offshore wind farm projects, the topic of atmospheric stability is unavoidable. In this post I’ll try to sum-up the way I usually handle it, since I am being asked often – I’d better have a text I can refer to :).

Sadly, I have waited a long time before looking into this seriously. That was certainly a mistake. I remember repeating the mantra that “strong winds are neutral” until at least year 2013 or so. Well, anyways, I have spent time catching up, and I am now routinely analysing atmospheric stability as part of my daily tasks.

Said shortly, assessing atmospheric stability means assessing how much atmospheric turbulence is created/dissipated by heat transfer, compared with mechanical friction. So, working with stability has a lot to do with trying to figure out which way (up or down) the heat flux goes, and what happens to the wind profile not only close to the surface but also higher-up. Luckily, we have now tons of publicly-available data to play with, I’ll show a few examples in this post.

If you are not familiar with the topic, you can find a great introduction in (Stull, 2017: Chap. 05 and Chap. 18). If you have heard of it already, you may know that the wind speed profile is typically well described, close to the surface and in flat terrain, by the Monin-Obukhov Similarity Theory (MOST) supplemented by some parametrisation of the roughness length. But somehow, even with practice, the topic is always a bit tricky, so in this post I have chosen to walk you through the way I understand it, step by step:

First, we’ll make a simple characterisation of the turbulent wind flow close to the surface, by splitting the flow velocity into a mean- and a fluctuating component, and derive the so-called turbulent stresses;
The next step will consist in making a simple model linking the vertical variation of the mean flow (the shear), to one of the turbulent stresses (the shear stress);
In the third step, we will add heat transfer to the model and discuss the MOST.
Finally, we will need to understand “how far up” our model works, and realise that it needs to be updated so as to reflect the real behavior of the wind field across the entire rotor (and not just close to the surface), in particular for stable conditions.

Step 1: Characterising turbulent stresses in the atmospheric boundary layer.

Let us look at wind measurements in the North Sea (horizontal wind speed at the IJmuiden met mast in the Dutch North Sea), and observe that the time series are turbulent, and that the level/kind of turbulence depend on the difference between the air- and the water temperature at the surface. When the air is warmer than the water, the level of turbulence relative to the mean speed is smaller than when the air is colder than the water.

Figure 1. Illustration of fluctuations (sampled at 4Hz) in the horizontal wind speed, over two distinct days at the top of the IJmuiden met mast (North Sea). The time series in blue comes from a summer day of 2015 where the air temperature was warmer than the sea surface temperature. The time series in red comes from a late autumn day of 2011 where the air temperature was colder than the sea surface temperature. When the water is warmer than the air, this adds energy to the flow in the form of heat, which is then transported vertically throughout the boundary layer, leading to larger and more powerful eddies than when the water is colder than the air. The time series in grey, visible on the 10-minute time series (rows 2 and 4) are from the other cup anemometer mounted at the top of the mast.

If we look at the variation of mean wind speed with elevation, it appears that when this temperature difference is small the mean wind speed profile is well represented by a logarithmic law. When the air is colder than the water the gradient $dU/dz$ (the shear) is smaller, and it is larger when the air is warmer than the water; also the corresponding profiles depart from a log law; see below.

Figure 2. Mean wind speed profile recorded at the Hollandse Kust Zuid A (HKZA) Floating LiDAR System (FLS) over the entire measurement campaigns (almost two years). The plot to the left shows the values with linear x-axis and logarithmic y-axis, while in the plot to the right both axes are in logarithmic scale. These subsets are defined according to the value of the difference between the air- and the water temperature.

So what is going on? If we consider a two-dimensional horizontal flow homogenous in the along-wind direction, and look at the subset of the data for which the temperature difference (and thereby, the heat flux) is small, the flow is characterised by:

the space- and time coordinates: vertical direction z and time t;
the air density ρ [kg/m³] and the air viscosity μ [kg/(m·s)]
the instantaneous horizontal wind speed u(t) [m/s];
the instantaneous vertical wind speed w(t) [m/s].

We are dealing with a turbulent flow, and it is practical to split the wind field in mean component (U,0) and a turbulent component (u‘,w’) and write:

u(t) = U + u‘(t)
w(t) = w‘(t)

We can use all of the above and write down the continuity- and conservation of momentum equation for the mean flow. This is often reffered to the as the Reynolds averaging of the Navier–Stokes equations, and leads to the emergence of a turbulent stresses (reynolds stresses):

In the along-wind direction: $\rho\overline{u^{'}u^{'}}$
In the vertical direction: $\tau = -\rho\overline{u^{'}w^{'}}$ .

The latter corresponds to the vertical transport of the turbulent fluctuations. The question is then: how to relate these microscale, turbulent properties of the flow, and the mean properties (such as the wind shear) ?

Step 2: Linking microscale turbulent stresses to mean wind shear

This is done by introducing a new parameter, called the mixing-length $l_{m}$ , which relates the fluctuation to the mean shear (in passing: the shear due to the flow velocity is neglected, due to the very small viscosity of the air):

$$u^\prime \sim l_{m} \frac{\partial U}{\partial z}$$

The mixing length is a characteristic distance over which the turbulent properties of a fluid particle remain constant. It is reasonable to assume that $u^{\prime} \sim w^\prime$ ; and thereby:

$$\overline{u^{\prime}w^{\prime}} \sim \left(l_{m}\frac{\partial U}{\partial z}\right)^{2}$$

For all practical purposes, $u_{*} = \sqrt{\overline{u^{\prime}w^{\prime}}}$ is called the “friction velocity“, and from experiments, close to the surface and when the atmosphere is well mixed:

$u_{*}$ is approximately constant;
the mixing length is proportional to the elevation $l_{m} \sim z.$

This results in the wind shear being inversely proportional to the elevation:

$$\frac{\partial U}{\partial z} \sim \frac{u_{*}}{z}$$

Experimentally, scientists and engineers have found that $l_{m} = \kappa z$ , where $\kappa$ is the Von Kármán constant; this results in:

$$\frac{\partial U}{\partial z} = \frac{u_{*}}{\kappa z}$$

The mean wind speed profile can then be derived by integrating the above expression, and introducing $z_{0}$ the elevation above the surface at which the speed reaches zero ( $z_{0}$ is called the “roughness length“):

$$U(z) = \frac{u_{*}}{\kappa}\ln\left(\frac{z}{z_{0}} \right)$$

At this stage, please note that, actually, this mixing-length assumption is not necessary for deriving the log law. This is explained in “The Logarithmic Wind Profile” (Tennekes, 1973), and nicely discussed in (Ortiz-Suslow et al., 2021).

It turns out that $z_{0}$ depends on $u_{*}$ : this is because the larger the waves get, the more “rough” the sea surface becomes. Also, for very small wind speeds, viscosity plays a role too; here is the form which such a relationship takes (from ECMWF):

$$z_{0} = \alpha_{M}\frac{\nu}{u_{*}} + \alpha_{ch}\frac{u_{*}^{2}}{g}$$

$\nu$ is the air viscosity;
$g$ is the acceleration of gravity;
$\alpha_{M}=0.11$ and $\alpha_{Ch}=0.018$ (Ch is for for Charnock, 1955) are constants derived from experiments.

So here we are, with a simple analytical model with which provide the wind speed at the elevation $z$ :

$$U(z) = \frac{u_{*}}{\kappa}\ln\left(\frac{z}{z_{0}(u_{*})} \right)$$

In this model, the wind speed is only a function of the elevation and the friction velocity. In turn $u_{*}$ and $z_{0}$ can be derived from a single wind speed value $U$ at elevation $z$ by solving the following implicit equation (easy to do with the fzero function in Matlab):

$$u_{*} = \frac{ U(z)\kappa} { \ln \left( \frac{z}{z_{0}(u_{*})} \right) }$$

How does this model compare with experiments? Well, but only for large wind speeds: see, the model is plotted in blue in the Figure below:

For small wind speeds, the measured shear is larger than the modelled shear;
For large wind speeds, the model fit well the measurements from 30 mMSL and up.
Note: the 4 mMSL measurements are recorded using a sonic anemometer mounted above the Floating LiDAR buoy, and it is possible that when the waves are large, the sensor sits “too low” in is somehow sheltered. It could also be that the sensor is wrong, but I favor the first explanation.

Figure 3. Comparison between measurements (black) and log-law model (blue), for the subset of the data where the difference between the air and the water is small.

So, why is that ? What makes the model fail for small wind speeds? As we shall see in the next Section, this is because until now we have neglected a key driver of atmospheric dynamics: buoyancy.

Step 3: Adding buoyancy – the Monin-Obukhov Similarity Theory

As we saw earlier in Figure 1 and Figure 2, turbulence and shear vary depending on the difference between the air- and the water temperature. Until now, in steps 1 and 2, we only considered mechanical turbulence created by the roughness of the surface, but in reality heat transfers between the fluid parcels lead to transfer of energy as well; this in turn influences the mixing-length and thereby the shear.

We see that when the water is warmer than the air, i.e. when there is an upward heat flux, the shear decreases. Coming back to the relationship between wind shear and momentum flux, this means that for such conditions the mixing length is larger:

$$\frac{\partial U}{\partial z} \sim \frac{u_{*}}{l_{m}} $$

And this make sense, thinking about what actually happens when the water is warmer than the air: this leads to convective (large) structures with length scales in the order of ~100 m (think about such convection cells for instance).

When warm air flows over colder water, there is no convection, and dominant microscale turbulent structures are smaller; i.e.: the mixing length is smaller and the shear becomes larger.

A major theoretical development on this topic came from (Richardson, 1920), who identified a criteria for the turbulence kinetic energy to increase in the atmosphere (the whole atmosphere, not just close to the surface), by comparing the mechanical shear-production/loss term of the turbulent kinetic energy budget equation and the buoyant production/consumption term; the ratio being the two terms has been coined the “Richardson number” $R_{i}$ . However, its formulation relies on the shear; sad (for us) because that is precisely what we are after. Yet, this formulation provides a first glimpse into the concept of atmospheric (static) stability: :

when $R_{i} <0$ the flow is unstable, and turbulent;
when $R_{i}$ is larger than a certain value (approx. 0.21), the flow stops being turbulent.

The second major theoretical development came from (Obukhov, 1946) and (Monin and Obukhov, 1954). In their work, the following two turbulent fluxes (microscale) are present: $\tau = \rho\overline{u^{\prime}w^{\prime}}$ which we saw earlier, and the buoyancy flux $q$ , which comes from the Reynolds averaging of the Navier Stokes equation when introducing buoyancy (variations of density due to potential temperature $\theta=\Theta+\theta^{\prime}$ , with $C_{p}$ the isobaric heat capacity for dry air):

$$q = \rho C_ {p} \overline{\theta^{\prime}w^{\prime}}$$

Monin and Obukhov then postulated that:

$\tau$ and $q$ are proportional respectively to the shear and the gradient of vertical potential temperature, via the mixing length (itself a function of $R_{i}$ );
close to the surface, the turbulent fluxes are constant.

They then postulated, based on dimensional arguments, that the gradients of wind speed and temperature could be expressed as functions of a new dimensionless parameter, $z/L$ , where $L$ is a characteristic length scaled defined as:

$$ L = – \frac{u_{*}^{3}T_{v}}{\kappa g \overline{\theta^{\prime}w^{\prime}}}$$

$L$ is the Obukhov length, $T_{v}$ is the mean temperature in the surface layer. For the wind speed gradient, the expression looks familiar:

$$\frac{\partial U}{\partial z} = \frac{u_{*}}{\kappa z}\phi_{M}\left(\frac{z}{L}\right)$$

$$U(z) = \frac{u_{*}}{\kappa} \left( \ln\left(\frac{z}{z_{0}(u_{*})}\right)+\psi_{M}\left(\frac{z}{L}\right)\right)$$

Yes! It’s basically the log-law in neutral conditions, modified by a factor which depends on the atmospheric stability. Scientists have spent a number of decades trying to measure and characterise experimentally the functions $\phi_{M}$ and $\psi_{M}$ (including folks from RISØ, Denmark). You can find expressions for this functions in various text books and articles. I use the ones in Section 3.3.2 of the ERA5 technical documentation part IV.

The ratio $z/L$ alone becomes a proxy for the atmospheric (dynamic) stability:

If $z/L$ is very small, then $\phi_M = 1$ ; the wind speed profile follows a log-law as we saw earlier; the effect of buoyancy can be neglected and the surface layer is reffered to as “neutral”.
If $L >0$ then the surface layer is stable ( $\phi_M >0$ ) and the shear is larger than in neutral conditions.
If $L <0$ then the surface layer is unstable ( $\phi_M >0$ ) and the shear is smaller than in neutral conditions.

Yayy, how nice. The only problem is that the computation of $L$ involves computing the covariance of high-frequency time series of $u$ , $w$ and $\theta$ : these are almost never available from measurements on commercial engineering projects. Also, models such as IFS or GFS do not resolve microscale turbulence, so how is it even possible to deal with the Monin-Obukhov similarity theory at a practical, engineering level?

So here is the secret:

from surface measurements, we derive a simplified version of the Richardson number called the bulk Richardson number $Ri_{b}$ ,
we use analytical formulas derived from experiments where both $Ri_{b}$ and $z/L$ have been measured conccurently to derive $z/L$ values.

This method has been used successfully in (Peña, 2008), where the author used the parametrisation provided in (Grachev and Fairall, 1997). The bulk Richardson number $Ri_{b}$ is defined as:

$$Ri_{b} = -\frac{gz\Delta\theta_{v}}{T(z)U(z)^{2}}$$

Where:

$g$ is the gravitational acceleration;
$z$ is an elevation above the sea surface (between 10 and 15 m);
$\Delta\theta_{v}$ is the absolute difference between the virtual potential temperature at elevation $z$ and the sea surface temperature (not the skin temperature, but the temperature slighlty below the surface);
$T$ and $U$ are the mean temperature and wind speeds.

Once we have obtained $Ri_{b}$ , $z/L$ is derived as follows:

$$ \frac{z}{L} = 10 Ri_{b} \text{ for unstable conditions}$$

$$ \frac{z}{L} = \frac{10 Ri_{b}}{1-5Ri_{b}} \text{ for stable conditions}$$

Using this method with the ERA5 data leads to the type of distributions below (with $z = 10 \text{mMSL}$ ), for a site in the North Sea. The classification is adapted from (Sathe, 2010).

Figure 4: Example of long-term distribution of hourly values of z/L computed from the ERA5 dataset using the method from (Peña, 2008).

Here are some simple, numerical examples with $z = 10 \text{m/s}$ , $\Delta\theta_{v} = 2 \text{K}$ and $T_{v} = 15 \text{K}$ :

With $U = 5 \text{m/s}$ , $Ri_{b} =0.0272$ and $L = 31.7 \text{m}$
With $U = 10 \text{m/s}$ , $Ri_{b} =0.0068$ and $L = 141.8 \text{m}$
With $U = 15 \text{m/s}$ , $Ri_{b} =0.0030$ and $L = 325.4 \text{m}$

Hey that’s my mantra: “strong winds are neutral“! Yes sure, strong winds are neutral, but weak winds matter too, in particular for doing energy yield assessments. Also, for wave-modelling purposes, or for analysing SAR wind data, it important to know about the stability for correctly computing the 10mMSL wind speed and the friction velocity. Let us take a 8 m/s wind speed at 10 m; then, solving the implicit equation mentioned earlier:

$u_{*} = 0.261 \text{m/s}$ for $L=50 \text{m}$
$u_{*} = 0.303 \text{m/s}$ for $L=-50 \text{m}$

The stress $\tau$ is proportional to the square of the friction velocity, so here the relative difference between the to stability conditions is about 35%!

Step 4: extending the model above the surface layer

Back to the MOST: this model works very well, but not for all elevations. As you can see in the Figure below, it predicts too large wind speeds above a certain elevation in stable conditions. This is because the assumptions behind the model break down:

the fluxes are not constant with elevation;
the mixing length stops increasing lineary, and after some elevation becomes constant.

Figure 5: This Figure shows, on the right, a comparison between the Monin-Obukhov Similarity Theory (MOST, dashed lines) and measurement data (dots). The full lines are mesoscale data from the Dutch European Wind Atlas. The color correspond to stability classes: unstable (red), neutral (black), stable (blue), very stable (purple), all (grey). The measurement location is the Hollandse Kust Zuid A floating LiDAR in the Dutch North Sea. The figure to the left shows the corresponding histograms of the height of the boundary layer, from ERA5 (where is is defined at the elevation where $Ri_{b}=0.25$ ; see Section 3.10.1 of ECMWF’s technical doc).

Figure 5: This Figure shows, on the right, a comparison between the Monin-Obukhov Similarity Theory (MOST, dashed lines) and measurement data (dots). The full lines are mesoscale data from the Dutch European Wind Atlas. The color correspond to stability classes: unstable (red), neutral (black), stable (blue), very stable (purple), all (grey). The measurement location is the Hollandse Kust Zuid A floating LiDAR in the Dutch North Sea. The figure to the left shows the corresponding histograms of the height of the boundary layer, from ERA5 (where is is defined at the elevation where $Ri_{b}=0.25$ ; see Section 3.10.1 of ECMWF’s technical doc).

These limitations are nicely explained in (Gryning, 2007): I have reproduced its Figure 1 below which depicts the evolution of the dominant length scale in the atmospheric boundary layer for neutral conditions. This papers explains very clearly how to extend the MOST above the surface layer, and all the way to the top of the boundary layer.

All of the details about can be found in the paper, and are very easy to implement. A comparison with measurement- and mesoscale data below shows that the model predicts very well the wind profile across a modern wind turbine rotor.

Figure 6: This Figure shows, on the right, a comparison between the Monin-Obukhov Similarity Theory extended above the surface layer in the model from Gryning (dashed lines) and measurement data (dots). The color correspond to stability classes: unstable (red), neutral (black), stable (blue), very stable (purple), all (grey). The full lines are mesoscale data from the Dutch European Wind Atlas. The measurement location is the Hollandse Kust Zuid A floating LiDAR in the Dutch North Sea. The figure to the left shows the corresponding histograms of the height of the boundary layer, from ERA5 (where is is defined at the elevation where $Ri_{b}=0.25$ ; see Section 3.10.1 of ECMWF’s technical doc).

We shall discuss in another post how this model can be used (spoiler: for deriving wind speed at any elevation, using CFSR/CFSv2 and ERA5 data). In the meantime, it is good to remember that sometimes reality is more complex than our simple models. For instance, warm air can be advected over water, causing stable conditions in altitude while near the surface winds remain neutral. This is discussed in the context of SAR measurements of wind farm wakes (loooooong wakes) in a previous post.

Comments/questions are welcome.

Rémi

Atmospheric stability is our friend (yes it is)

Step 1: Characterising turbulent stresses in the atmospheric boundary layer.

Step 2: Linking microscale turbulent stresses to mean wind shear

Step 3: Adding buoyancy – the Monin-Obukhov Similarity Theory

Step 4: extending the model above the surface layer

Published by Rémi

1 Comment