Transport Error#

A footprint says where a measurement’s air came from, according to one meteorological analysis. The analysis has wind errors, so the modelled enhancement has an error too. PYSTILT estimates it with the method of Lin and Gerbig (2005). You run the particles a second time with a random wind error added, one that has the statistics of the meteorology’s errors. The perturbed particles spread further. The extra spread of the modelled enhancement across the particles is the transport-error variance. An inversion uses it as the transport part of each observation’s model-data mismatch.

Run with wind errors#

The perturbed run is a variant of its own. It uses the same receptors and meteorology, with four wind-error settings added. The names are the same as in STILT-R.

variants:
  hrrr: {}                 # the unperturbed run
  hrrr-err:                # the same run with perturbed winds
    siguverr: 2.6          # wind error standard deviation, m/s
    tluverr: 260           # its correlation time, min
    zcoruverr: 450         # its vertical correlation length, m
    horcoruverr: 14        # its horizontal correlation length, km
    grid: null             # particles only

Every receptor then has a second simulation, model.simulations[rid, "hrrr-err"], whose particles saw the perturbed winds. You can add the error variant to a finished project, and only the new simulations run. The error run costs as much as the original. grid: null skips the footprint, which the error calculation does not need. Keep a grid if you want to plot a perturbed footprint.

Choosing the settings#

The correlation scales matter as much as the standard deviation. HYSPLIT decorrelates the wind error over time (tluverr) and over the distance a particle travels (horcoruverr). If the scales are too short, the error turns into noise that averages out along each trajectory. The perturbed particles then spread no further than the unperturbed ones, and there is nothing to measure. At 10 m/s, for example, a particle crosses a 5 km correlation length in eight minutes.

Derive the four values from wind observations for your meteorology and region, as Lin and Gerbig did (Wind Error Statistics). Do not copy them from another analysis, and do not pick small scales to be safe. The values above are for HRRR over the Salt Lake Valley.

HYSPLIT also has mixed-layer height errors (sigzierr, tlzierr, horcorzierr), set the same way. They have little effect. HYSPLIT multiplies each particle’s footprint by its own random factor and leaves the path alone. Because the factors are independent between particles, they average out over a few thousand particles. For a mixed-layer error that every particle shares, use ziscale instead (Mixed-layer height).

Several realizations#

One error run is one random draw of the wind error. To average over several draws, set realizations on the error variant:

hrrr-err:
  siguverr: 2.6
  # ...
  realizations: 4          # hrrr-err-0 .. hrrr-err-3

The variant then runs four times per receptor, as the simulations hrrr-err-0 to hrrr-err-3. The unperturbed hrrr run is not repeated. A variant with realizations always gets numbered names, even with realizations: 1. Each realization needs a different random draw, and there are two ways to get one:

  • krand: 4, the default. HYSPLIT seeds each run from the clock. The realizations differ, but a rerun gives different ones.

  • krand: 2 with a seed. Realization k runs with seed + k, and a rerun reproduces every realization exactly.

PYSTILT rejects any other setting when it reads the config, because every realization would get the same draw. Each realization is its own simulation, so a rerun only does the missing ones. You can raise realizations later and only the new realizations run.

Pass all the realizations to transport_error() as a list. It averages them before taking the difference:

sims = model.simulations.sel(receptor=rid)
err = transport_error(
    sims[rid, "hrrr"].trajectories.data,
    [t.data for t in sims.sel(variant="hrrr-err").trajectories.load().values()],
    flux,
)
err.realizations  # 4

More realizations help less than you might expect. They reduce the noise of the perturbed runs, but every realization is compared with the same unperturbed run, so the noise of the estimate falls by at most a factor of √2 (see Notes). Realizations are worth their cost when a single run’s variance is within a factor of two of its noise. When variance is far below noise, check the wind-error scales instead.

The modelled enhancement#

The enhancement a footprint predicts for a surface flux field is the footprint times the flux, summed over the grid. flux is an xarray.DataArray on lat and lon, with an optional time dimension. It is sampled at the footprint’s cell centres, and footprint cells outside the flux field add nothing. If the flux cells are much smaller than the footprint cells, regrid the flux first.

flux = xr.open_dataarray("ch4_flux.nc")          # µmol m⁻² s⁻¹ on lat/lon
foot = model.simulations[rid, "hrrr"].footprint
enhancement = foot.enhancement(flux)             # ppm per footprint time step
total = float(enhancement.sum())

A flux in µmol m⁻² s⁻¹ times a footprint in ppm per (µmol m⁻² s⁻¹) gives ppm.

The emission error#

The flux field’s own uncertainty gives another error on the enhancement. Take sigma, a field of one standard deviation per cell in the flux’s units. Put it on the footprint grid and compute two limits:

import numpy as np

f = foot.integrate_over_time()
s = sigma.reindex(lat=f["lat"], lon=f["lon"], method="nearest")
err_correlated = float((f * s).sum())                   # every cell errs the same way
err_independent = float(np.sqrt(((f * s) ** 2).sum()))  # each cell on its own

The first assumes the errors in all cells move together, so it is an upper bound. X-STILT reports this one as the emission error (cal.emiss.err), with sigma from the spread of several inventories. The second assumes the cells are independent, which gives a lower bound. Real inventory errors are correlated over some distance, so the true value lies in between. Computing it needs the covariance of the flux errors, which is the prior error covariance of an inversion. If you set up the inversion with fips, InverseProblem.prior_obs_error is that covariance carried to every observation (the footprint matrix times the prior covariance times its transpose). The square root of its diagonal is each observation’s emission error. Add it in quadrature to the transport and retrieval errors to get the total error of a modelled value.

The transport error#

transport_error() takes a receptor’s unperturbed and perturbed particle tables and the flux field:

import pandas as pd
from stilt.observations import transport_error

rows = []
for sim in model.simulations.sel(variant="hrrr"):
    err = model.simulations[sim.id.receptor, "hrrr-err"]
    result = transport_error(
        sim.trajectories.data,
        err.trajectories.data,
        flux,
        transforms=sim.config.transforms,
        context=sim.transform_context(),
    )
    rows.append({"receptor": sim.receptor.id, "enhancement": result.enhancement,
                 "variance": result.variance, "noise": result.noise})
errors = pd.DataFrame(rows).set_index("receptor")

result.variance is the transport-error variance of the modelled enhancement, in the enhancement’s units squared. result.sd is its square root. Pass the footprint’s transforms and the simulation’s context so the particles are weighted the way the footprint weights them. That covers the averaging kernel (including one from a per-receptor table), pressure weighting, and lifetime decay. For a tower receptor there is nothing to pass.

Is the estimate meaningful?#

variance is the difference of two sample variances, each from a few thousand particles. It is noisy and is sometimes negative. Two checks tell you whether a value means anything:

  • result.noise is the standard deviation variance would have with no wind error at all. A variance within two or three times noise is not resolved.

  • Over many receptors, take the median of the signed variance, by hour, season, or site. Do not clip each value at zero first, because that turns noise into a positive error. result.sd treats a negative variance as zero, so use it only once the variance is resolved.

The signal is largest when turbulence spreads the particles least, at night and in winter. On a convective afternoon the unperturbed particles already fill the boundary layer. The wind error adds little, and the estimate is usually within its noise.

What the numbers mean#

result.levels shows the calculation for each release level. A point receptor has one level.

  • height is the level’s mean release height in metres, and n its number of particles.

  • weight is the level’s share of the particles.

  • mean_orig and mean_err are the mean enhancement per particle without and with the perturbation. Their weighted sums are result.enhancement and result.enhancement_perturbed.

  • var_orig and var_err are the variances of the per-particle enhancement, and dvar is their difference. This is Lin and Gerbig’s equation 4 for the level.

  • sd_trans is the square root of dvar, keeping its sign.

result.enhancement is close to what the footprint gives. The two differ a little because the footprint smooths the particles onto its grid.

A column receptor’s particles are grouped into levels release-height bins of equal width, 20 by default. A multipoint or slant receptor with no more than levels release heights gets one level per height. The levels are then combined with an exponential error correlation in the vertical:

\[\sigma^2 = \sum_i w_i^2\, \mathrm{dvar}_i + \sum_{i \ne j} w_i w_j s_i s_j\, e^{-|h_i - h_j| / L}\]

Here w is weight, s is sd_trans, h is height, and L is length_scale, 356 m by default (X-STILT’s value). With length_scale=None the levels are independent.

percentile=0.99 reproduces X-STILT’s trimming (Wu et al., 2018). It drops the top 1 % of particles in each level before taking the variance, which damps the few particles that pass over a point source at the cost of a small bias. X-STILT also fits a line through the levels where dvar is positive. PYSTILT does not, because fitting only those levels biases the result upward: under pure sampling noise it reports a positive error at every level.

The method measures how much the wind error moves particles between flux cells. It says little when the flux field is uniform, and it is only as good as the wind statistics you gave the run. It covers transport only. Emission and retrieval errors are separate terms. Wind errors also move the trajectory endpoints, and with them the background. Pass a background field to include that part (Background).

Mixed-layer height#

A real error in the mixed-layer height is shared. Every particle in the valley sees the same layer, too shallow or too deep. ziscale models this. It multiplies HYSPLIT’s mixed-layer height by one factor for every particle. Runs with factors above and below 1.0 show how sensitive the enhancement is to the mixed layer. Declare them as variants of the same project:

variants:
  hrrr: {}
  hrrr-zi06: {ziscale: 0.6}   # every hour of the run; a list gives one factor per hour
  hrrr-zi14: {ziscale: 1.4}

Before you run them:

  • Set krand: 2 and a seed in the defaults. Every variant then draws the same turbulence, and the difference between variants is the mixed layer alone. Without a seed each run adds its own sampling noise.

  • HYSPLIT applies the minimum mixing depth kmix0 (150 m by default) after the factor. A mixed layer already at that floor is not lowered further.

  • HYSPLIT holds at most 150 hourly factors. A single ziscale value is repeated for every hour, so it needs abs(n_hours) <= 150. For a longer run, give a list. Hours past the end of the list are unscaled.

How much it matters depends on the receptor. A column spans the mixed layer, so a change in its depth mostly moves influence around inside the column, and the column enhancement changes little. A surface receptor has no such averaging and can change by tens of percent. Test it for your own receptors.

The spread between these runs is a sensitivity. To turn it into an error you need to know how far the meteorology’s mixed-layer height is from the real one. One way is to compare it with radiosonde profiles, analysed with the bulk Richardson method HYSPLIT uses by default (kmixd: 3).

Notes#

  • Lin and Gerbig derived their scales from variograms of analysed minus radiosonde winds, and got about 120 km, 4 hours, and 900 m for an 80 km analysis. For a 3 km model such as HRRR the horizontal scale is much shorter (Wind Error Statistics).

  • With krand: 4 the clock gives only about 5000 distinct seeds, so two realizations can be identical. With a few realizations this is very unlikely. With a hundred it is more likely than not.

  • With krand: 2, realization 0 uses seed itself, as STILT-R’s error run does.

  • noise comes from splitting the unperturbed particles into random halves and treating one half as the perturbed run (noise_splits times, 16 by default). With N realizations the noise is sqrt((1 + 1/N) / 2) times that of a single run, which never falls below 1/sqrt(2). result.noise already includes this factor.