Use power low instead of log law for extrapolating wind speed

[lukas-rokka]

Use power low instead of log law for extrapolating wind speed to hub height.

Description

The extrapolate_wind_speed function uses the power law and the grid averages forecast surface roughness (fsr): wnd_spd = ds[from_name] * ( np.log(to_height / ds["roughness"]) / np.log(from_height / ds["roughness"]) ) Using the ERA5 fsr with the log law results in too large correction, especially in the Nordics where the fsr is often high.
At height above the blending height it is better to use the power law. That looks like follow: wnd_spd = ds["wnd100m"] * (to_height / 100)**a. The wind sheer exponent a can be derived from the wind speed at 10m and 100m: a = np.log(ds["wnd10m"] / ds["wnd100m"]) / np.log(10 / 100). Tested it and it give same results as doing actual vertical interpolations with MetPy using the wind components from the pressure level data (how the ERA5 wind at 100 meters is orginally calculated as well).

For turbines hub heights lower than the blending height it could be useful to use the log law to account for site specific roughness length. Below the equations from the IFS CY41R2 docs used to derive the 10m wind speed in ERA5 (always calculated assuming fsr of 0.03 to better represent conditions at weather stations)

The IFS docs mention two blending heights, 40 and 75 meters (I guess it actually might be dynamic and site specific). But a procedure I use my self is to first use the above-mentioned power law to calculate a wind speed at a blending height of 70 meters, then the log law using a site specific fsr (for now, just some blunt estimates based on the ERA5 grid average fsr)

[euronion]

Thanks @lukas-rokka .

Could you share some of the results you obtained with this method compared to the current method or MetPy? Also if you have more information on the method that would be welcome :).

Do you think this should be the one-and-only method, or could/should this be an alternative to the current approach?

For reference IFS CY41R2.

[lukas-rokka]

In this post it says that the wind at 100 meters is calculated using interpolation on the log scale https://confluence.ecmwf.int/pages/viewpage.action?pageId=155343870. Doing linear interpolation on the log scale is mathematically equivalent to the power law mentioned above.

This paper uses the power law and the wind speed at 10 and 100 meters to calculate the exponent: https://doi.org/10.3390/en14144169. Guess there will be a slight error introduced when extrapolating above 100 meters as the ERA5 10m is calculated using the log law and a fsr of 0.03, but shouldn’t be too large.

Guess it could be good having both extrapolation methods available. This site says that the log law is appropriate for the lowest 100 meters. So, this issue should maybe be labelled feature and not as a bug 😊

Haven't done any well worked through testing, but here a code that can be used to calculate the wind speed at certain heights from pressure levels. Doesn't return exactly the same result as ERA5 100m (maybe because ERA5 seem to use quadratic interpolation accoding to the docs)

import xarray as xr
import metpy.calc as mpcalc
import numpy as np
ds = xr.open_datset( #a dataset containing the geopotential and the u and v wind components at 1000, 975 and 950 pressure levels. 
ds["altitude"] = mpcalc.geopotential_to_height(ds.z)
ds["altitude"].values = ds["altitude"].values # remove unit

lat =59; lon=20 
height_to = 100 #+ 116 # hub_height + geopotential height at surface level

a1000_975 = np.log((ds.ws.sel(level=1000))/ds.ws.sel(level=975)) / np.log((ds.altitude.sel(level=1000))/ds.altitude.sel(level=975))
a975_950 = np.log((ds.ws.sel(level=975))/ds.ws.sel(level=950)) / np.log((ds.altitude.sel(level=975))/ds.altitude.sel(level=950))

ws100a = ds.ws.sel(level=975) * ((height_to)/ds.altitude.sel(level=975))**a1000_975
ws100b = ds.ws.sel(level=950) * ((height_to)/ds.altitude.sel(level=950))**a975_950
ws100 = np.where(ds.altitude.sel(longitude=lon, latitude=lat, level=1000) > 0, ws100a.sel(longitude=lon, latitude=lat), ws100b.sel(longitude=lon, latitude=lat))

MetPy log_interpolate_1d can also be used https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.log_interpolate_1d.html.

[euronion]

Thanks! Feature it is then. We could also have a parameter interpolation_method which takes ['automatic', 'log', 'power'] as arguments and if automatic is selected, then the appropriate logic based on the hub_height is selected.

We'd be happy to receive a PR for this :) Else we'll probably have a look at it in May. Next weeks will be too busy with other things.

[lukas-rokka]

Task 2.3 in the new PECDv4 suggest using the power law as well https://climate.copernicus.eu/sites/default/files/2022-02/C3S2_412_Volume_II_final.pdf. It also says that there might be alpha coefficient (aka wind shear exponents) made available in CDS. So might be better to wait for what is being delivered from that project.

I’m not that familiar with either PR or the structure of atlite so not sure I’m right person for the task. But basically, the steps would be to

1. download also 10 m wind speeds
2. calculate time varying alpha coefficient from wind 10m and 100m. ds["a"] = np.log(ds["wnd10m"] / ds["wnd100m"]) / np.log(10 / 100). (wind speed at 10 m can be discarded after this step).
3. An extrapolate_wind_speed_power_law function, something like wnd_hub = ds["wnd100m"] * (to_height / 100)**ds["a"]
4. Maybe a wrapper function for doing automatic, log law or power law extrapolation.

Additionally or alternatively use alpha coefficient provide by PECDv4.