Wind Error Statistics#

A transport-error run (Transport Error) needs four numbers that describe the meteorology’s wind errors. They are the standard deviation of the error and how long, how high, and how far it stays correlated. Lin and Gerbig (2005) derive them by comparing the analysed wind with observed winds. This page does the same for your meteorology and region. Derive your own values instead of copying them from another analysis or place. The correlation scales decide whether the perturbation has any effect at all.

What the numbers are#

For each wind component, the error at an observation is the analysed wind minus the observed wind. Its standard deviation is siguverr.

The correlation scales come from the variogram. At a separation h, the variogram is half the mean squared difference between the errors of pairs of observations h apart. For an error with an exponential correlation it is

\[\gamma(h) = \sigma^2 \left(1 - e^{-h/l}\right)\]

where σ is the standard deviation and l the correlation scale. Fitting this curve with σ fixed gives l. Pairs of heights give zcoruverr, pairs of times give tluverr, and pairs of stations give horcoruverr.

Each scale needs data that resolves it. Radiosonde profiles resolve height, but launches are twelve hours apart, so they cannot measure a time scale of a few hours. Hourly surface stations resolve time and horizontal distance, but they only sample the bottom of the layer. The recipe below takes the standard deviation and the vertical scale from radiosondes, and the time and horizontal scales from surface stations.

Step 1: the errors#

Sample the meteorology at the observations with arlmet. Winds in ARL files on projected grids are stored relative to the grid, so ask for earth-relative components (earth_relative=True, arlmet 0.1.0a8 or later). Radiosondes report geopotential height, so sample in metres above sea level:

import arlmet
import pandas as pd

sondes = pd.read_csv("slc_2024_sondes.csv")   # time, lon, lat, height (m MSL), elevation, u, v
points = sondes.rename(columns={"height": "z"})
met = arlmet.sample_points(hrrr_files, points, ["UWND", "VWND"],
                           z_kind="msl", earth_relative=True)
upper = pd.DataFrame({
    "time": sondes["time"],
    "height": sondes["height"] - sondes["elevation"],   # m above ground
    "u_err": met["UWND"] - sondes["u"],
    "v_err": met["VWND"] - sondes["v"],
})

Surface stations report wind at a height above ground, 10 m for a standard anemometer. Sample the 10 m wind, or the lowest model level:

# stations: time, site, lon, lat, u, v
points = stations.assign(z=10.0)
met = arlmet.sample_points(hrrr_files, points, ["U10M", "V10M"],
                           z_kind="agl", earth_relative=True)
surface = pd.DataFrame({
    "time": stations["time"], "site": stations["site"],
    "lon": stations["lon"], "lat": stations["lat"],
    "u_err": met["U10M"] - stations["u"],
    "v_err": met["V10M"] - stations["v"],
})

Reading the observations is up to you. The Salt Lake Valley numbers below used the sonde reader in lair.

For radiosondes anywhere in the world, the usual source is NOAA’s Integrated Global Radiosonde Archive (IGRA2). siphon reads it (pip install siphon). IGRA2 names a WMO station by its two-letter country code, M, and the WMO number padded to eight digits. Salt Lake City, WMO 72572, is USM00072572. The full list is igra2-station-list.txt in the NCEI archive. This builds the sondes table above:

import datetime as dt
from siphon.simplewebservice.igra2 import IGRAUpperAir

levels, launches = IGRAUpperAir.request_data(
    [dt.datetime(2024, 1, 1), dt.datetime(2024, 12, 31, 23)], "USM00072572"
)
launches = launches.drop_duplicates("date")[["date", "latitude", "longitude"]]
levels = levels.merge(launches, on="date")
ground = levels["lvltyp2"] == 1                  # IGRA2's surface level
elevation = levels["height"].where(ground).groupby(levels["date"]).transform("max")
sondes = pd.DataFrame({
    "time": levels["date"],
    "lon": levels["longitude"], "lat": levels["latitude"],
    "height": levels["height"],                  # geopotential height, m
    "elevation": elevation,
    "u": levels["u_wind"], "v": levels["v_wind"],
}).dropna()

Each call downloads the station’s whole record and keeps the dates you ask for, so ask for the full period once, not one launch at a time. time is the nominal launch hour, 00 or 12 UTC. The balloon goes up as much as an hour earlier. The last line drops levels without a height, which IGRA2 has on some significant levels. It also drops launches with no surface level, since their heights above ground are unknown.

Step 2: the scales#

variogram() builds the empirical variogram. It takes the errors, the coordinate the separation is measured in, and a group label. Only points with the same label are paired. The coordinate can also be two columns of longitude and latitude, and the separation is then the great-circle distance in kilometres. fit_variogram() fits the exponential model with σ fixed.

import numpy as np
from stilt.observations import fit_variogram, variogram

def scale(errors, lag, group, bins):
    table = variogram(errors, lag, group=group, bins=bins)
    return fit_variogram(table["lag"], table["gamma"], sigma=errors.std()).length

layer = upper[upper["height"].between(0, 3000)]        # m above ground
minutes = (surface["time"] - pd.Timestamp("2000-01-01")) / pd.Timedelta("1min")
components = ("u_err", "v_err")

siguverr = np.mean([layer[c].std() for c in components])
zcoruverr = np.mean([scale(layer[c], layer["height"], layer["time"],
                           range(0, 3100, 100)) for c in components])
tluverr = np.mean([scale(surface[c], minutes, surface["site"],
                         range(0, 14401, 60)) for c in components])
horcoruverr = np.mean([scale(surface[c], surface[["lon", "lat"]], surface["time"],
                             range(0, 51)) for c in components])

The recipe makes these choices, and you can change any of them:

layer is the height range whose errors matter for your transport. Use the boundary layer and a little above for a surface receptor, and more for a column.
The vertical variogram pairs levels within one launch (group is the launch time), in 100 m bins up to 3 km.
The time variogram pairs hours at one station, out to ten days.
The horizontal variogram pairs stations at one time, out to 50 km.
Each fit fixes σ at the standard deviation of the errors it was built from.
u and v are done separately and averaged.
The mean error (the bias) is ignored, as in Lin and Gerbig. The perturbation does not represent it.

The four values go into config.yaml under the same names, or into ErrorParams. Look at the fits before you trust them:

import matplotlib.pyplot as plt

table = variogram(layer["u_err"], layer["height"], group=layer["time"],
                  bins=range(0, 3100, 100))
fit = fit_variogram(table["lag"], table["gamma"], sigma=layer["u_err"].std())
plt.scatter(table["lag"], table["gamma"], s=8)
plt.plot(table["lag"], fit(table["lag"]))

Seasonal and diurnal values are worth a look. Filter both tables and run the recipe again. For the Salt Lake Valley the horizontal scale ranged from 9 km in summer to 20 km in winter.

Without surface stations, the time scale can only come from pairs of launches at the same height. That tells you whether the errors are still correlated twelve hours later, and little else. horcoruverr then has to come from another source.

Salt Lake Valley, HRRR, 2024#

Radiosondes from Salt Lake City and 20 surface stations in the valley, compared with the 3 km HRRR analysis, gave:

Component	σ [m/s]	l_z [m]	l_x [km]
u	2.37	365	12.1
v	2.87	548	16.5

These became the values in the Transport Error example: siguverr: 2.6, zcoruverr: 450, horcoruverr: 14 and tluverr: 260. The hourly stations gave time scales of 207 min for u and 221 min for v, and the sondes, extrapolated, 212 to 312 min. The recipe above gives 2.5 m/s, 480 m, 214 min and 14 km. For comparison, Lin and Gerbig found about 120 km, 4 hours and 900 m for an 80 km analysis over the eastern United States.