Wind profile program

Logarithmic wind profile

Theory

Atmospheric bound layer books, such as `The atmospheric boundary layer' by J.R. Garratt, will show that during the day and near the surface the velocity is expected to vary according to a logarithmic equation

  u = A log ( (z + h) / z0)

where z is height from surface, z₀ is the roughness length, h is a displacement height and A = u_* / κ - where u_* is the friction velocity and &kappa is the Von Karman constant.

This equation is derived because a convective boundary layer is expected during the day, and this implies constant stress which leads to the equation above. There is a displacement height, h, because there may be some roughness elements, such as waves or bushes, that mean that the variation of stress with height does not settle to a constant until somewhere above the surface.

If there's some stability, there will be Monin-Obukhov correction to this, but this will generally be small and the equation above should be a better fit than a simple linear interpolation. It is nighttime when the boundary layer is expected to be most stable, although I would still expect the equation above to be a better fit than a simple linear interpolation. Under almost any circumstances the wind speed will increase most rapidly at the surface and this will be better represented by a log equation than a linear interpolation.

Fitting a log equation to the wind speed

The equation above has three unknowns, A, h and z₀, which can be solved because we have three wind speeds at 2m, 10m and 50m, which we'll call speed2m, speed10m and speed50m. It can be shown that

  A = speed50m / log((50.0 + h) / z0) 

and that

  z0 = Exp ( ( speed50m * log (2.0 + h) + speed2m * log (50.0 + h) ) / 
      (50.0 - 2.0) ).

Hence, solving h will give us A and z₀. Using the equations at all three heights, can produce an implicit equation for h

  (speed10m - speed2m) * log ( (50.0 + h)/(2.0 + h) ) - 
      (speed50m - speed2m) * log ( (10.0 + h)/(2.0 + h) = 0.

Let the function f(h) be given by

  f(h) = (speed10m - speed2m) * log ( (50.0 + h)/(2.0 + h) ) - 
          (speed50m - speed2m) * log ( (10.0 + h)/(2.0 + h) = 0.

and this implies that at the correct displacement height that f(h) = 0. And the Newton-Raphson method can be used to iteratively find the solution to h,

  hn+1 = hn - f(hn) / f'(hn)

around the neighbourhood of the solution. It can be shown that the gradient of f(h) is given by

  f'(h) = ( speed50m * (50.0 + h) * (10.0 - 2.0) - 
            speed10m * (10.0 + h) * (50.0 - 2.0) +
            speed2m  * (2.0  + h) * (50.0 - 10.0) ) /
            ( (50.0 + h) * (10.0 + h) * (2.0 + h) )

which is defined provided that h > -2.0. And this is positive provided that h < h_max, where

  hmax = ( (speed50m - speed2m) * (10.0 - 2.0) * 50.0 -
                      (speed10m - speed2m) * (50.0 - 2.0) * 10.0 ) /
                    ( (speed10m - speed2m) * (50.0 - 2.0) -
                      (speed50m - speed2m) * (10.0 - 2.0) )

so h needs to be kept between -2.0 and h_max. Hence the iterative equation to solve h becomes

 hn+1 = hn - (50.0 + h) * (10.0 + h) * (2.0 + h) *
                   [ (speed10m - speed2m) * log ( (50.0 + h)/(2.0 + h) ) - 
                     (speed50m - speed2m) * log ( (10.0 + h)/(2.0 + h) ) ] /
                    [ speed50m * (50.0 + h) * (10.0 - 2.0) - 
                      speed10m * (10.0 + h) * (50.0 - 2.0) +
                      speed2m  * (2.0  + h) * (50.0 - 10.0) ].

This is iterated until |h_n+1 - h_n| < 0.0001 or until 80 loops is reached.

Failure to solve this occurs when speed10m ~ speed50m, or when speed2m ~ speed10m with a large speed50m.

At the Reading location (-0.9685° East, 51.445° North) there were 13 profiles where a log equation was not found when using all the data we had (at 6:30 on 21 May 1983, 10:30 on 31 May 1983, 16:30 on 23 August 1983, 11:30 on 27 August 1983, 11:30 27 August 1984, 7:30 on 10 September 1988, 15:30 on 12 February 2004, 7:30 on 14 October 2004, 15:30 on 28 March 2005, 23:30 on 2 June 2008, 10:30 on 31 August 2008, 20:30 on 8 August 2009 and 18:30 on 21 June 2010).