library(ggplot2)
te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
ink = "#16241d", paper = "#f5f4ee")
theme_te <- function() {
theme_minimal(base_size = 12) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "#e7e6dc"),
plot.background = element_rect(fill = "#f5f4ee", colour = NA),
panel.background = element_rect(fill = "#f5f4ee", colour = NA),
plot.title = element_text(face = "bold", colour = te_pal$ink),
axis.title = element_text(colour = "#2c3a31"),
axis.text = element_text(colour = "#2c3a31"))
}Night-time flux and the u-star threshold
R
eddy covariance
carbon flux
data quality
ecology tutorial
A simulated year of half-hourly eddy covariance data measures the friction velocity filter: what it recovers, what it discards and what it costs a budget.
Corrected 4 November 2026: the post said that Aubinet (2008) concluded advection cannot be reliably measured with the instrumentation most towers carry; his abstract supports only that direct advection measurements are highly site dependent and that the friction velocity threshold alone may be insufficient, and the text now says that. It also called the simulated true exchange a sink of moderate size for a temperate peatland, which nothing in the post supports; it now says only that the sink is set by the simulation’s parameters. The introduction said the post runs the estimator practitioners use; it now says a moving point estimator of the kind Papale et al (2006) describe, since this post runs the test once on the whole year.
Corrected 14 October 2026: the post put the colder soil of calm nights down to clear skies, called its whole-year moving point test the form standardised by Papale et al and their highest-quarter rule a convention of other implementations, and said that Papale et al resample whole nights and that Massman and Lee find the corrections interact; it now puts the colder soil mostly down to the season, says that Papale et al run the test on each quarter of the year and filter at the highest of the four thresholds, where this post runs it once on the whole year, that they resample single half hours, and that Massman and Lee note that many of the shortcomings of eddy covariance coincide at night. The over-imputation of the single mean fill is now an average over the discarded half hours rather than the same on every one of them.
Corrected 1 October 2026: the post said that noise only ever ends the moving point walk too soon, put the nonzero threshold on the unattenuated flux down to residual temperature gradients and called the whole gap at the estimated threshold systematic; it now says that noise can also let the walk run on, that the rule can only return the mean of one of its bins, and that part of that gap is the luck of this noise draw. The first figure’s caption and alt text now say that the attenuated points tend to sit low within their nights, where they had said nothing in their position marked them out, the temperature gradient across the pooled bins is put at more than a doubling of respiration rather than tens of per cent, a share of night-time half hours is no longer called a share of nights, and the dropped values are described as positive on average rather than all positive.
A flux tower stands on a temperate fen at forty-seven degrees north, sonic anemometer and gas analyser on a short mast a couple of metres above a sedge canopy. It has been running for a year. The logger writes one record every half hour, so the year is seventeen thousand five hundred and twenty rows of net ecosystem exchange, air and soil temperature, incoming radiation and friction velocity. The question the tower was put up to answer is whether the fen took up carbon over the year and how much.
The sign convention first, because it causes more confusion than anything else here. Net ecosystem exchange is written from the atmosphere’s point of view: a negative flux is carbon leaving the air and entering the ecosystem, a positive flux is carbon going the other way. So daytime half hours in the growing season are negative, every night-time half hour is positive, and an annual sum that comes out negative describes a sink.
Eddy covariance works by measuring the covariance between vertical wind speed and gas concentration, which means it needs vertical wind. On a calm clear night there is very little. The air near the ground stratifies, turbulence collapses, and the carbon dioxide respired by the soil and the plants does not reach the sensor: it pools in the canopy air below it, or drains sideways down whatever slope exists, and is measured hours later or not at all. Goulden et al (1996) set this out in the first careful accuracy audit of a long eddy covariance record, and the standard response has been the same ever since: measure the friction velocity, decide on a threshold below which the night-time flux is not to be trusted, and throw those half hours away. This post simulates a year in which the true flux is known at every half hour, applies that attenuation, and runs a moving point estimator of the kind Papale et al (2006) describe, to see how close it gets, what it costs in data and what it does to the annual number.
What this post is and is not next to
The nearest neighbour on this site is the litter decomposition cluster, and the difference is worth stating because the two look superficially similar and fail in completely different ways. Mass loss and the carbon budget also converts a measurement into a carbon flux, but there the measurement is a bag on a balance: a discrete object, weighed on eight occasions over five years, with corrections that are multiplicative and systematic, to do with what the mass is made of and what fell into the bag. Here nothing is weighed and the year is an integral of seventeen thousand half hours. The failure mode is not a wrong conversion factor applied to a correct number; it is a subset of the record, selected by weather rather than at random, in which the instrument was reading low, plus the decision about which subset that was.
This is the first of four posts on eddy covariance. This one stops at the threshold: it estimates it, prices it and leaves gaps in the record. Filling those gaps is gap filling a flux time series; splitting the filled record into gross primary production and ecosystem respiration is partitioning net flux into GPP and respiration; auditing the finished figure is checking an annual flux budget.
A year of half hours
The drivers come first: solar geometry from latitude and day of year, cloudiness as a day-to-day autocorrelated series, air temperature with a seasonal cycle, a diurnal cycle and synoptic wobble, and a soil temperature that is a damped and lagged version of it. One feature of the temperature field is put in deliberately and matters later. Daytime heating is stronger under clear skies and night-time cooling is too: on a clear night the surface loses heat to space, the loss accumulates, and the near-surface air and the top of the soil keep falling until dawn. That accumulation is written as a term in the square root of the hours since sunset, the usual shape for radiative cooling under a stable layer.
n_half <- 17520
tix <- seq_len(n_half)
doy <- (tix - 1) %/% 48 + 1
hod <- ((tix - 1) %% 48) * 0.5 + 0.25
lat_deg <- 47
declin <- 23.44 * sin(2 * pi * (284 + doy) / 365) * pi / 180
sin_elev <- sin(lat_deg * pi / 180) * sin(declin) +
cos(lat_deg * pi / 180) * cos(declin) * cos((hod - 12) * 15 * pi / 180)
is_night <- sin_elev <= 0
rl <- rle(is_night)
hrs_dark <- 0.5 * sequence(rl$lengths) * rep(as.integer(rl$values), rl$lengths)
night_id <- (tix + 23) %/% 48 + 1
set.seed(20260802)
cloud <- plogis(as.numeric(stats::filter(rnorm(365, 0, 0.8), 0.62,
method = "recursive")) - 0.15)
clear <- 1 - rep(cloud, each = 48)
ppfd <- pmax(0, 2100 * sin_elev) * (0.22 + 0.72 * clear)
seas_t <- 9 + 9.5 * sin(2 * pi * (doy - 112) / 365)
syn <- rep(as.numeric(stats::filter(rnorm(365, 0, 1.7), 0.72,
method = "recursive")), each = 48)
anom <- (2 + 7 * clear) * pmax(0, sin_elev) - 2.1 * clear * sqrt(hrs_dark)
tair <- seas_t + syn + anom
ts_slow <- as.numeric(stats::filter(0.004 * (seas_t + syn), 0.996,
method = "recursive",
init = seas_t[1] + syn[1]))
tsoil <- ts_slow + 0.42 * anom
print(round(c(half_hours = n_half, nights = max(night_id),
night_fraction = mean(is_night),
longest_night_hours = max(hrs_dark)), 3)) half_hours nights night_fraction longest_night_hours
17520.0 366.0 0.5 16.0 print(round(c(air_min = min(tair), air_max = max(tair),
soil_min = min(tsoil), soil_max = max(tsoil),
peak_light = max(ppfd)), 2)) air_min air_max soil_min soil_max peak_light
-8.94 28.50 -4.42 22.40 1592.51 The record covers 17520 half hours over 366 night windows, the first and last of them partial, and 50 per cent of it falls with the sun below the horizon, the longest night running to 16 hours. Air temperature spans -8.9 to 28.5 degrees C and soil temperature -4.4 to 22.4, with peak photosynthetic photon flux density of 1593 umol/m2/s.
The true fluxes are built from two standard response functions. Ecosystem respiration follows the temperature response of Lloyd and Taylor (1994), driven by soil temperature, and gross primary production follows a rectangular hyperbola in light whose asymptote is scaled by a seasonal greenness curve. Net ecosystem exchange is respiration minus production, so it is positive at night and negative in the middle of a summer day.
r_ref <- 2.2
e0_lt <- 300
reco <- r_ref * exp(e0_lt * (1 / (283.15 - 227.13) -
1 / ((tsoil + 273.15) - 227.13)))
green <- 1 / (1 + exp(-(doy - 108) / 11)) * 1 / (1 + exp((doy - 288) / 14))
green <- 0.02 + 0.98 * green
a_max <- 18 * green
alph <- 0.038
gpp <- a_max * alph * ppfd / (a_max + alph * ppfd)
gpp[ppfd <= 0] <- 0
nee_true <- reco - gpp
conv <- 1800 * 12.011e-6
c_reco <- sum(reco) * conv
c_gpp <- sum(gpp) * conv
c_true <- sum(nee_true) * conv
print(round(c(reco_ref_umol = r_ref, activation_e0 = e0_lt,
quantum_yield = alph, amax_peak = max(a_max)), 4))reco_ref_umol activation_e0 quantum_yield amax_peak
2.2000 300.0000 0.0380 17.9739 print(round(c(annual_reco = c_reco, annual_gpp = c_gpp,
annual_nee = c_true,
night_mean_flux = mean(nee_true[is_night]),
day_min_flux = min(nee_true)), 3)) annual_reco annual_gpp annual_nee night_mean_flux day_min_flux
865.119 1096.581 -231.462 1.817 -9.273 Over the year the fen respires 865 gC/m2 and fixes 1097 gC/m2, so the true net exchange is -231.46 gC/m2/yr, a net sink set by the parameters of this simulation. Night-time flux averages 1.817 umol/m2/s and the strongest midsummer uptake reaches -9.27 umol/m2/s. That annual figure is the thing every method below is trying to recover, and in the field it does not exist.
The defect: what a calm night does to the measurement
Friction velocity is the square root of the momentum flux and it is the field measure of how much mechanical turbulence there is. The series here has four parts: a daytime convective boost that follows solar elevation, a synoptic wind component that persists for a few days at a time, a seasonal roughness term tied to the greenness curve, since a sedge canopy in July is rougher than the same fen in February, and a stability penalty on clear nights that deepens as the night goes on. The last two are what link turbulence to temperature. A clear calm night is both the least turbulent and the coldest; a windy cloudy one is neither. That link is physics rather than an artefact of the simulation, and it is why the estimator in the next section is more complicated than it first looks as though it needs to be.
set.seed(20260803)
w_syn <- rep(as.numeric(stats::filter(rnorm(365, 0, 0.32), 0.70,
method = "recursive")), each = 48)
eps_u <- as.numeric(stats::filter(rnorm(n_half, 0, 0.11), 0.86,
method = "recursive"))
ustar <- exp(log(0.40) + 0.55 * pmax(0, sin_elev)^0.6 + 0.32 * green + w_syn -
0.55 * clear * is_night - 0.055 * clear * hrs_dark + eps_u)
u_mid <- 0.16
u_w <- 0.040
u_true <- u_mid + u_w * log(99)
att <- ifelse(is_night, plogis((ustar - u_mid) / u_w), 1)
nee_att <- nee_true * att
print(round(quantile(ustar[is_night], c(0.05, 0.25, 0.5, 0.75, 0.95)), 3)) 5% 25% 50% 75% 95%
0.116 0.214 0.311 0.440 0.782 print(round(c(ramp_midpoint = u_mid, ramp_width = u_w,
true_threshold = u_true,
transfer_at_midpoint = plogis(0),
transfer_at_threshold = plogis((u_true - u_mid) / u_w)), 4)) ramp_midpoint ramp_width true_threshold
0.1600 0.0400 0.3438
transfer_at_midpoint transfer_at_threshold
0.5000 0.9900 print(round(c(mean_night_transfer = mean(att[is_night]),
mean_flux_lost = mean(nee_true[is_night] - nee_att[is_night]),
night_below_threshold = mean(ustar[is_night] < u_true),
carbon_lost = sum(nee_true * (1 - att)) * conv), 4)) mean_night_transfer mean_flux_lost night_below_threshold
0.8473 0.2177 0.5811
carbon_lost
41.2389 The attenuation is a logistic ramp in friction velocity rather than a step. Below a midpoint of 0.16 m/s more than half the respired carbon fails to reach the sensor; well above it the flux passes intact. Defining the true threshold as the friction velocity at which 99 per cent of the flux gets through puts it at 0.3438 m/s, and that is the number every estimate below is compared against. A step would make the estimation problem trivial and would also be wrong, since the collapse of turbulent transport is gradual.
Over the year 58.1 per cent of night-time half hours sit below that threshold, the average night-time half hour delivers 84.7 per cent of its true flux, and the carbon that never reaches the sensor totals 41.2 gC/m2 across the year. Measurement noise goes on top. Real flux noise is not constant: its scale grows with the size of the flux, as Richardson et al (2006) showed across a range of towers, so the standard deviation here is an offset plus a fraction of the flux magnitude.
set.seed(20260804)
sig <- 0.45 + 0.12 * abs(nee_att)
nee_obs <- nee_att + rnorm(n_half, 0, sig)
nu <- ustar[is_night]
nt <- tsoil[is_night]
nf <- nee_obs[is_night]
nid <- night_id[is_night]
print(round(c(noise_floor = 0.45, noise_slope = 0.12,
mean_sigma = mean(sig), sigma_night = mean(sig[is_night]),
sigma_day = mean(sig[!is_night])), 4))noise_floor noise_slope mean_sigma sigma_night sigma_day
0.4500 0.1200 0.7506 0.6419 0.8592 print(round(c(sd_night_flux = sd(nf), mean_night_flux = mean(nf),
annual_sum_noise_sd = sqrt(sum(sig^2)) * conv), 4)) sd_night_flux mean_night_flux annual_sum_noise_sd
1.4728 1.5934 2.3037 The noise standard deviation on a night-time half hour averages 0.642 umol/m2/s, against a night-time flux that averages 1.593 with a standard deviation of 1.473. Summed over the whole year the noise contributes a standard deviation of only 2.304 gC/m2/yr to the annual total, because seventeen thousand independent errors mostly cancel. The section on how well the threshold is known measures the threshold uncertainty against the same quantity, estimated there by redrawing the noise and passing it through the gap filling.
The fall-off is not visible in the record
The first thing to establish is that nothing in the time series announces the problem. Twelve days of summer half hours, with the badly attenuated ones marked, look like twelve days of ordinary flux data.
Two measurements say the same thing with numbers: the size of the loss against the size of the ordinary variation, and how much of the variance in night-time flux friction velocity accounts for, before and after soil temperature has been taken out of it.
loss_night <- mean(nee_true[is_night] - nee_att[is_night])
fit_t <- lm(nf ~ nt)
fit_u <- lm(nf ~ nu)
fit_p <- lm(residuals(fit_t) ~ nu)
nm_f <- as.numeric(tapply(nee_obs[is_night], nid, mean))
nm_u <- as.numeric(tapply(ustar[is_night], nid, mean))
calm <- nm_u < u_true
gap_nightly <- mean(nm_f[!calm]) - mean(nm_f[calm])
nm_t <- as.numeric(tapply(tsoil[is_night], nid, mean))
nm_cool <- as.numeric(tapply((tsoil - ts_slow)[is_night], nid, mean))
nm_doy <- as.numeric(tapply(doy[is_night], nid, max))
cold_q <- nm_doy <= 90 | nm_doy > 273
soil_gap <- mean(nm_t[!calm]) - mean(nm_t[calm])
cool_gap <- mean(nm_cool[!calm]) - mean(nm_cool[calm])
print(round(c(mean_flux_lost = loss_night, sd_night_flux = sd(nf),
loss_over_sd = loss_night / sd(nf)), 4))mean_flux_lost sd_night_flux loss_over_sd
0.2177 1.4728 0.1478 print(round(c(r2_soil_temperature = summary(fit_t)$r.squared,
r2_ustar = summary(fit_u)$r.squared,
r2_ustar_after_temperature = summary(fit_p)$r.squared), 4)) r2_soil_temperature r2_ustar
0.6732 0.1484
r2_ustar_after_temperature
0.0755 print(round(c(calm_nights = sum(calm), other_nights = sum(!calm),
mean_flux_calm = mean(nm_f[calm]),
mean_flux_other = mean(nm_f[!calm]),
apparent_gap = gap_nightly,
actual_loss = loss_night,
gap_over_loss = gap_nightly / loss_night), 4)) calm_nights other_nights mean_flux_calm mean_flux_other apparent_gap
187.0000 179.0000 1.2131 2.4618 1.2487
actual_loss gap_over_loss
0.2177 5.7347 print(round(c(soil_colder_on_calm = soil_gap,
from_clear_sky_cooling = cool_gap,
calm_in_oct_to_mar = mean(cold_q[calm]),
other_in_oct_to_mar = mean(cold_q[!calm])), 4)) soil_colder_on_calm from_clear_sky_cooling calm_in_oct_to_mar
4.4286 0.2629 0.6578
other_in_oct_to_mar
0.3352 The attenuation removes 0.2177 umol/m2/s from the average night-time half hour, against a night-time standard deviation of 1.473: a ratio of 0.148. Soil temperature accounts for 67.3 per cent of the variance in night-time flux and friction velocity for 14.8 per cent on its own, falling to 7.5 per cent once temperature has been removed first. Anyone looking at the record without a specific reason to bin it by turbulence will see temperature and nothing else.
The nightly aggregate makes the trap sharper. Averaging each night and splitting the 366 nights at the true threshold gives 187 calm nights averaging 1.213 umol/m2/s and 179 others averaging 2.462, an apparent gap of 1.249. That gap is 5.73 times the flux the attenuation actually removed. Most of it is temperature, and little of that comes from clear skies: the soil is 4.43 degrees C colder on calm nights than on the others, only 0.26 of it from the extra cooling of clear nights, and the rest is mostly the season, with 66 per cent of the calm nights falling between October and March against 34 per cent of the others. A reader who took that gap for the size of the problem would be out by nearly a factor of six, in a direction that is not obvious.
The signature, and where it hides
The way to see the attenuation is to sort the night-time half hours by friction velocity, cut them into equal-sized bins and take the mean flux in each. If turbulence is doing nothing to the measurement, the bin means are flat. If it is, they climb from the left and then level off, and the levelling-off point is the threshold. Aubinet et al (2000) made that plot part of the standard EUROFLUX processing chain and it has been the diagnostic ever since.
bin_table <- function(us, fx, tp, n_bin = 20) {
bn <- cut(us, quantile(us, seq(0, 1, length.out = n_bin + 1)),
include.lowest = TRUE, labels = FALSE)
data.frame(bin = seq_len(n_bin),
ustar = as.numeric(tapply(us, bn, mean)),
flux = as.numeric(tapply(fx, bn, mean)),
soil = as.numeric(tapply(tp, bn, mean)),
n = as.numeric(table(bn)))
}
pool_tab <- bin_table(nu, nf, nt)
t_class <- cut(nt, quantile(nt, seq(0, 1, length.out = 7)),
include.lowest = TRUE, labels = FALSE)
in_cls <- which(t_class == 4)
cls_tab <- bin_table(nu[in_cls], nf[in_cls], nt[in_cls])
print(round(pool_tab[c(1, 4, 8, 12, 16, 20), ], 3)) bin ustar flux soil n
1 1 0.096 0.272 5.490 438
4 4 0.183 0.908 4.185 438
8 8 0.264 1.451 5.395 438
12 12 0.341 1.864 6.845 438
16 16 0.462 2.061 7.705 438
20 20 1.025 2.828 12.160 438print(round(cls_tab[c(1, 4, 8, 12, 16, 20), ], 3)) bin ustar flux soil n
1 1 0.073 0.214 7.842 73
4 4 0.114 0.511 8.082 73
8 8 0.189 1.131 7.842 73
12 12 0.282 1.715 7.635 73
16 16 0.408 1.795 8.135 73
20 20 0.952 1.962 9.026 73print(round(c(pooled_soil_range = diff(range(pool_tab$soil)),
class_soil_range = diff(range(cls_tab$soil)),
pooled_flux_low = pool_tab$flux[1],
pooled_flux_high = mean(pool_tab$flux[16:20]),
class_flux_low = cls_tab$flux[1],
class_flux_high = mean(cls_tab$flux[16:20])), 3))pooled_soil_range class_soil_range pooled_flux_low pooled_flux_high
8.608 1.422 0.272 2.379
class_flux_low class_flux_high
0.214 1.867 Pooling every night-time half hour in the year gives a curve that rises across the whole range, from 0.272 umol/m2/s in the lowest bin to 2.379 in the top five, and never plateaus. The soil temperature column says why: mean soil temperature in the pooled bins spans 8.61 degrees C, because the seasonal roughness of the canopy has loaded the high friction velocity bins with summer nights, so the curve is measuring the season as much as the turbulence. Restricting to one of six soil temperature classes shrinks that spread to 1.42 degrees C and the shape changes: the flux climbs from 0.214 umol/m2/s to a plateau at 1.867 and stops. That is the signature, and temperature classing is what makes it appear.
The moving point test
The estimator in general use is the moving point test of Gu et al (2005), in the form standardised by Papale et al (2006). Split the night-time record into temperature classes. Inside each class, sort by friction velocity and cut into twenty equal-sized bins. Walk up the bins and take the first one whose mean flux reaches ninety-nine per cent of the mean of all the bins above it: that bin’s friction velocity is the class threshold. Discard any class in which friction velocity and temperature are correlated beyond a set limit, on the grounds that the classing has failed there. Take the median of what is left. Papale et al (2006) do all of this separately for each quarter of the year and filter the full record at the highest of the four quarterly thresholds; the code below runs it once on the whole year, a simplification the section on the honest limit comes back to.
mpt_bin <- function(us, fx, n_bin = 20, crit = 0.99) {
if (length(us) < 5 * n_bin) return(NA_real_)
bn <- cut(us, quantile(us, seq(0, 1, length.out = n_bin + 1)),
include.lowest = TRUE, labels = FALSE)
m_f <- as.numeric(tapply(fx, bn, mean))
m_u <- as.numeric(tapply(us, bn, mean))
for (i in seq_len(n_bin - 1)) {
if (m_f[i] >= crit * mean(m_f[(i + 1):n_bin])) return(m_u[i])
}
m_u[n_bin]
}
mpt_classed <- function(us, fx, tp, n_class = 6, n_bin = 20,
crit = 0.99, r_max = 0.4) {
cl <- cut(tp, quantile(tp, seq(0, 1, length.out = n_class + 1)),
include.lowest = TRUE, labels = FALSE)
th <- rep(NA_real_, n_class)
rho <- rep(NA_real_, n_class)
for (k in seq_len(n_class)) {
i <- which(cl == k)
rho[k] <- cor(us[i], tp[i])
if (abs(rho[k]) <= r_max) th[k] <- mpt_bin(us[i], fx[i], n_bin, crit)
}
list(th = th, rho = rho, med = median(th, na.rm = TRUE))
}
r_max_used <- 0.4
cc <- mpt_classed(nu, nf, nt, r_max = r_max_used)
u_hat <- cc$med
print(round(rbind(threshold = cc$th, ustar_temp_cor = cc$rho), 4)) [,1] [,2] [,3] [,4] [,5] [,6]
threshold 0.2047 0.2419 0.2501 0.9524 0.2732 0.3250
ustar_temp_cor 0.0166 0.2140 -0.1736 0.1603 -0.3244 0.0631print(round(c(estimated_threshold = u_hat, true_threshold = u_true,
error = u_hat - u_true,
error_pct = 100 * (u_hat / u_true - 1),
classes_used = sum(!is.na(cc$th))), 4))estimated_threshold true_threshold error error_pct
0.2617 0.3438 -0.0821 -23.8909
classes_used
6.0000 The six classes return 0.205, 0.242, 0.25, 0.952, 0.273 and 0.325 m/s. One of them is absurd: the fourth class returns 0.952, which is above almost every night-time friction velocity in the record and means the rule never triggered inside that class at all. This is why the summary across classes is a median rather than a mean, and it is the single most useful piece of engineering in the whole procedure. The median comes out at 0.2617 m/s against a true 0.3438, an error of -0.0821 m/s or -23.9 per cent.
The estimate is low, and three comparisons separate the candidate reasons: what the rule would converge on with infinite clean data, what it returns on the attenuated flux with the noise switched off, and what it returns on the true unattenuated flux.
u_grid <- seq(0.05, 0.7, by = 0.001)
plateau <- vapply(u_grid, function(x) mean(att[is_night][nu > x]), numeric(1))
u_target <- u_grid[which(plogis((u_grid - u_mid) / u_w) >= 0.99 * plateau)[1]]
u_clean <- mpt_classed(nu, nee_att[is_night], nt)$med
u_notrb <- mpt_classed(nu, nee_true[is_night], nt)$med
print(round(c(true_threshold = u_true,
population_target_of_rule = u_target,
rule_on_noise_free_attenuated = u_clean,
rule_on_noisy_attenuated = u_hat,
rule_on_unattenuated_flux = u_notrb), 4)) true_threshold population_target_of_rule
0.3438 0.3370
rule_on_noise_free_attenuated rule_on_noisy_attenuated
0.3522 0.2617
rule_on_unattenuated_flux
0.1352 print(round(c(criterion_effect = u_target - u_true,
binning_effect = u_clean - u_target,
noise_effect = u_hat - u_clean), 4))criterion_effect binning_effect noise_effect
-0.0068 0.0152 -0.0905 The rule’s own target is not the true threshold. Because the ninety-nine per cent criterion is measured against the mean of the higher bins, and those bins are themselves attenuated by a hair, the value the rule converges on is 0.337 m/s rather than 0.3438, an offset of -0.0068. Applied to the attenuated flux with no measurement noise the classed rule returns 0.3522, which overshoots that target by 0.0152 and is close enough to the truth to be called correct.
The noise is what breaks it. Adding the measurement error moves the estimate from 0.3522 to 0.2617, a shift of -0.0905 m/s, which is the whole of the error. The mechanism is in the shape of the rule: it takes the first bin that crosses the criterion, so a bin whose mean is pushed up by chance stops the walk early, twenty bins each get a chance to end the search, and on balance the search ends too soon, although a bin pushed down by chance can also let the walk run on, as it did in the fourth class above. A first-crossing rule is biased downwards by noise, and more data fixes it only slowly, because bin means sharpen with the square root of the record length while the net error stays one-sided. That is one of the arguments Barr et al (2013) make for replacing the walk with a fitted change point, which uses all the bins at once. Run the same classed rule on the true unattenuated flux, where there is no threshold to find, and it returns 0.1352 m/s: not zero, because the rule can only return the mean friction velocity of one of its bins, and the two classes that set the median both stopped at their lowest bin, the earliest the walk can stop. That is close to the lowest value the procedure can return, so it is not inventing thresholds out of nothing.
What the temperature classing buys
Skipping the classing is the obvious simplification, and it is a disaster here. Running the same walk on the pooled night-time record gives an answer, and the answer is meaningless.
u_pool <- mpt_bin(nu, nf)
u_pool_true <- mpt_bin(nu, nee_true[is_night])
cor_all <- cor(nu, nt)
print(round(c(pooled_estimate = u_pool,
classed_estimate = u_hat,
true_threshold = u_true,
pooled_bias = u_pool - u_true,
pooled_over_classed = u_pool / u_hat), 4)) pooled_estimate classed_estimate true_threshold pooled_bias
1.0253 0.2617 0.3438 0.6815
pooled_over_classed
3.9184 print(round(c(cor_ustar_soil_all_night = cor_all,
max_abs_within_class = max(abs(cc$rho)),
mean_abs_within_class = mean(abs(cc$rho)),
night_data_above_pooled = mean(nu >= u_pool)), 4))cor_ustar_soil_all_night max_abs_within_class mean_abs_within_class
0.2781 0.3244 0.1587
night_data_above_pooled
0.0169 print(round(c(pooled_on_unattenuated = u_pool_true,
top_bin_ustar = pool_tab$ustar[20]), 4))pooled_on_unattenuated top_bin_ustar
1.0253 1.0253 The pooled walk returns 1.0253 m/s, which is not an overestimate of the threshold but the rule failing to trigger at all and falling out of the loop at the top bin, hence a value equal to that bin’s mean friction velocity of 1.0253. Applying it would keep 1.69 per cent of the night-time record. Run it on the true unattenuated flux and it returns 1.0253 as well: the pooled curve gives the same answer when there is nothing to find.
The correlation between friction velocity and soil temperature across all night-time half hours is 0.278. Inside the six temperature classes the largest absolute correlation is 0.324 and the average 0.159, all below the limit of 0.4 that Papale et al (2006) set for discarding a class. The classing has not abolished the confounding; it has cut it to a size the rule can survive. A correlation of 0.278 sounds mild, and the reason something that mild does so much damage is the size of the criterion: the rule is testing for a one per cent difference in mean flux, while the temperature gradient across the pooled bins is 8.61 degrees C, worth more than a doubling of respiration at this temperature sensitivity. A confounder does not have to be strong. It only has to be stronger than one per cent.
What the filter costs in data
Two fractions are worth reporting and they are routinely confused. The filter removes night-time half hours only, so its effect on the night-time record is roughly twice its effect on the record as a whole.
gapi <- is_night & ustar < u_hat
keep_n <- is_night & !gapi
frac_all <- mean(gapi)
frac_night <- mean(nu < u_hat)
frac_at_true <- mean(nu < u_true)
run_len <- rle(gapi)
long_runs <- run_len$lengths[run_len$values]
print(round(c(discarded_half_hours = sum(gapi),
fraction_of_all = frac_all,
fraction_of_night = frac_night,
fraction_of_night_at_true_threshold = frac_at_true,
nights_wholly_lost = sum(tapply(gapi[is_night], nid, all))), 4)) discarded_half_hours fraction_of_all
3239.0000 0.1849
fraction_of_night fraction_of_night_at_true_threshold
0.3697 0.5811
nights_wholly_lost
40.0000 print(round(c(gap_runs = length(long_runs),
median_run_half_hours = median(long_runs),
longest_run_half_hours = max(long_runs),
mean_true_flux_discarded = mean(nee_true[gapi]),
mean_true_flux_retained = mean(nee_true[keep_n])), 4)) gap_runs median_run_half_hours longest_run_half_hours
412.0000 4.0000 32.0000
mean_true_flux_discarded mean_true_flux_retained
1.4480 2.0334 At the estimated threshold the filter throws away 3239 half hours, which is 18.49 per cent of the year and 36.97 per cent of the night-time record. Had it been set at the true threshold it would have taken 58.11 per cent of the night-time half hours instead, so the estimator’s downward error is also, quietly, a decision to keep an extra 21.13 per cent of night-time data that should have gone.
The gaps are not scattered. They come in 412 runs with a median length of 4 half hours and a longest run of 32, and 40 nights lose every one of their half hours: a gap of one half hour can be interpolated, a whole week of calm November nights cannot. The true flux on the discarded half hours averages 1.448 umol/m2/s and on the retained ones 2.033, so what the filter removes is not a random sample of the night but systematically the colder, quieter part of it, with genuinely lower respiration. Any method that fills those gaps with an average of the surviving data will put back too much carbon.
The annual number
Now the sum. Five versions, all of the same year.
c_nofilt <- sum(nee_obs) * conv
c_drop <- sum(nee_obs[!gapi]) * conv
imp_flat <- mean(nee_obs[keep_n])
c_flat <- (sum(nee_obs[!gapi]) + sum(gapi) * imp_flat) * conv
br_t <- quantile(tsoil[is_night], seq(0, 1, length.out = 11))
br_t[1] <- -Inf
br_t[length(br_t)] <- Inf
cl_t <- cut(tsoil, br_t, labels = FALSE)
annual_cls <- function(thr, obs = nee_obs) {
gi <- is_night & ustar < thr
nk <- is_night & !gi
mm <- as.numeric(tapply(obs[nk], factor(cl_t[nk], levels = seq_len(10)), mean))
if (anyNA(mm)) mm[is.na(mm)] <- mean(obs[nk])
(sum(obs[!gi]) + sum(mm[cl_t[gi]])) * conv
}
c_cls <- annual_cls(u_hat)
budget <- data.frame(
method = c("known truth", "no filter at all", "filter, gaps dropped",
"filter, one mean for all gaps", "filter, mean within temperature class"),
gC = c(c_true, c_nofilt, c_drop, c_flat, c_cls))
budget$error <- budget$gC - c_true
print(round(budget[, c("gC", "error")], 2)) gC error
1 -231.46 0.00
2 -273.08 -41.61
3 -336.37 -104.91
4 -196.46 35.00
5 -235.38 -3.92print(budget$method)[1] "known truth"
[2] "no filter at all"
[3] "filter, gaps dropped"
[4] "filter, one mean for all gaps"
[5] "filter, mean within temperature class"print(round(c(imputed_value = imp_flat,
true_mean_on_gaps = mean(nee_true[gapi]),
over_imputation = imp_flat - mean(nee_true[gapi])), 4)) imputed_value true_mean_on_gaps over_imputation
1.9979 1.4480 0.5499 Against a true -231.46 gC/m2/yr, the unfiltered sum comes out at -273.08, overstating the sink by 41.61 gC/m2/yr or 18 per cent. That is the attenuation arriving in the answer: carbon that left the fen but was never measured leaving it.
Filtering and then dropping the discarded half hours from the sum gives -336.37, which is worse than doing nothing at all by a wide margin, for a reason that is arithmetic rather than ecological. A sum over a subset is not an estimate of the sum over the whole; dropping a value from a total is imputing zero for it, and the values being dropped are night-time releases, positive on average. The filter identified the bad half hours correctly and then the sum threw away the carbon along with them.
Replacing each gap with the mean of the retained night-time data gives -196.46, an error of 35 gC/m2/yr in the other direction. This is the honest name for what dropping was trying to be: mean imputation. It is better than zero imputation by a factor of 3, and it is still wrong, for the reason the previous section flagged. The imputed value is 1.998 umol/m2/s while the true mean on the gaps is 1.448, an over-imputation averaging 0.55 umol/m2/s over all 3239 half hours.
The last row is the crudest possible improvement: sort the retained night-time data into ten soil temperature classes and fill each gap with the mean of its own class. That gives -235.38, an error of -3.92 gC/m2/yr, so a single covariate used in the bluntest possible way removes almost the whole of the imputation error. There is no comparison of filling methods here, and that is the next post’s job; the point is that the residual error after filtering is an imputation problem rather than a filtering problem, and that the choice of fill changes the annual number by more than the filter does.
How well is the threshold known
The threshold is a statistic and it has a sampling distribution. The resampling unit has to be the night, not the half hour: half hours inside one night share a weather state, a temperature and a turbulence regime, and resampling them individually would treat a year as if it held seventeen thousand independent observations of turbulence rather than a few hundred. Papale et al (2006) resample single half hours instead, and report that of all the processing choices in a flux chain the friction velocity correction contributes most to the uncertainty of an annual figure.
set.seed(20260805)
nights <- unique(nid)
n_night <- length(nights)
by_night <- split(seq_along(nid), nid)
n_boot <- 300
boot_u <- numeric(n_boot)
for (b in seq_len(n_boot)) {
pick <- unlist(by_night[as.character(sample(nights, n_night, replace = TRUE))],
use.names = FALSE)
boot_u[b] <- mpt_classed(nu[pick], nf[pick], nt[pick])$med
}
annual_flat <- function(thr, obs = nee_obs) {
gi <- is_night & ustar < thr
nk <- is_night & !gi
(sum(obs[!gi]) + sum(gi) * mean(obs[nk])) * conv
}
boot_flat <- vapply(boot_u, annual_flat, numeric(1))
boot_cls <- vapply(boot_u, annual_cls, numeric(1))
print(round(c(replicates = n_boot, nights_resampled = n_night,
point_estimate = u_hat, boot_mean = mean(boot_u),
boot_sd = sd(boot_u),
boot_q05 = unname(quantile(boot_u, 0.05)),
boot_q95 = unname(quantile(boot_u, 0.95)),
true_threshold = u_true), 4)) replicates nights_resampled point_estimate boot_mean
300.0000 366.0000 0.2617 0.2776
boot_sd boot_q05 boot_q95 true_threshold
0.0199 0.2463 0.3076 0.3438 print(round(c(flat_sd = sd(boot_flat), flat_q05 = unname(quantile(boot_flat, 0.05)),
flat_q95 = unname(quantile(boot_flat, 0.95)),
class_sd = sd(boot_cls), class_q05 = unname(quantile(boot_cls, 0.05)),
class_q95 = unname(quantile(boot_cls, 0.95))), 3)) flat_sd flat_q05 flat_q95 class_sd class_q05 class_q95
9.045 -204.782 -176.314 1.537 -237.267 -232.328 Over 300 resamples of the 366 nights the threshold has a standard deviation of 0.0199 m/s and a fifth to ninety-fifth percentile range of 0.2463 to 0.3076. The true value, 0.3438, lies above the ninety-fifth percentile. That is the first thing to say about the bootstrap: it is a sampling interval around a biased point estimate, and it does not contain the answer.
Pushing each of those 300 thresholds through the annual sum gives a standard deviation of 9.04 gC/m2/yr with the flat imputation and 1.54 with the temperature-class imputation. The threshold’s contribution to the budget’s uncertainty is not a property of the threshold. It is a property of the threshold and the gap filling together, and it differs between the two by a factor of 5.9.
Which matters more to the annual figure, the threshold or the random measurement noise? Both are measurable here: redraw the noise on a fixed threshold, and compare the spread of the annual sum with the spread the bootstrap produced.
set.seed(20260806)
n_noise <- 200
noise_flat <- numeric(n_noise)
noise_cls <- numeric(n_noise)
for (b in seq_len(n_noise)) {
obs_b <- nee_att + rnorm(n_half, 0, sig)
noise_flat[b] <- annual_flat(u_hat, obs_b)
noise_cls[b] <- annual_cls(u_hat, obs_b)
}
set.seed(20260807)
n_both <- 150
both_flat <- numeric(n_both)
both_u <- numeric(n_both)
for (b in seq_len(n_both)) {
obs_b <- nee_att + rnorm(n_half, 0, sig)
both_u[b] <- mpt_classed(nu, obs_b[is_night], nt)$med
both_flat[b] <- annual_flat(both_u[b], obs_b)
}
sys_gap <- annual_flat(u_hat) - annual_flat(u_true)
print(round(c(noise_only_flat = sd(noise_flat), noise_only_class = sd(noise_cls),
threshold_only_flat = sd(boot_flat),
threshold_only_class = sd(boot_cls)), 3)) noise_only_flat noise_only_class threshold_only_flat
2.653 2.600 9.045
threshold_only_class
1.537 print(round(c(ratio_flat = sd(boot_flat) / sd(noise_flat),
ratio_class = sd(boot_cls) / sd(noise_cls),
noise_through_threshold = sd(both_flat),
threshold_sd_from_noise = sd(both_u)), 4)) ratio_flat ratio_class noise_through_threshold
3.4088 0.5910 9.2540
threshold_sd_from_noise
0.0171 print(round(c(annual_at_estimate = annual_flat(u_hat),
annual_at_truth = annual_flat(u_true),
systematic_gap = sys_gap,
gap_over_boot_sd = abs(sys_gap) / sd(boot_flat)), 3))annual_at_estimate annual_at_truth systematic_gap gap_over_boot_sd
-196.462 -161.736 -34.727 3.839 Under the flat imputation it is the threshold, by a clear margin: a standard deviation of 9.04 gC/m2/yr against 2.65 from the measurement noise, a ratio of 3.41. Under the temperature-class imputation it is the noise: 1.54 against 2.6, a ratio of 0.59. The familiar claim that the threshold dominates the random error is true of this simulation under a crude gap fill and false under a slightly less crude one, so stating it without saying which gap fill is meant says nothing.
Two further numbers matter more than the ratio. Redrawing the noise and re-estimating the threshold from the redrawn data gives a spread of 9.25 gC/m2/yr, with the threshold itself varying by 0.0171 m/s across those draws. So the measurement noise is not a separate small term to be added to the threshold uncertainty: a large part of the threshold uncertainty is measurement noise, arriving through the estimator rather than through the sum.
And the error from using the estimated threshold is larger than both. Evaluating the annual sum at the estimated threshold gives -196.46 gC/m2/yr and at the true threshold -161.74, a gap of -34.73 gC/m2/yr, which is 3.8 times the bootstrap standard deviation. Not all of that gap is systematic: this draw’s threshold is among the lowest of the 150 noise redraws above, which average 0.294 m/s, and at that average threshold the gap is -20.97 gC/m2/yr, still 2.3 times the bootstrap standard deviation. A published uncertainty built from resampling alone would quote a range that excludes the value the correct threshold produces. Note also which of the two is closer to the truth: the biased threshold, at -196.46 against a true -231.46, beats the correct one at -161.74, because a filter that discards more data gives the biased imputation more to be wrong about. Two errors partly cancelling is not a method.
The honest limit
Everything above rests on one property of the simulation that no tower has: the attenuation is a known function of a single variable. Friction velocity was the only thing that decided how much flux got through, so a filter on friction velocity could in principle be exactly right, and the only question was whether the estimator found the number.
A real calm night is not like that. The carbon that fails to reach the sensor has at least three destinations and the filter addresses one of them. Some is stored in the canopy air below the instrument and released in a burst at sunrise, which a storage term from a concentration profile can partly recover and which single-point storage estimates handle badly. Some drains downslope as a density current, and on anything other than flat ground that advective loss is not a function of friction velocity at all: Aubinet (2008) analyses the nocturnal problem, finds direct advection measurements highly site dependent, and suggests that the friction velocity threshold alone may be insufficient. Massman and Lee (2002) set the friction velocity filter alongside the other corrections a flux record needs and note that many of the shortcomings of eddy covariance coincide at night. Some of it simply goes missing.
The seasonality in this simulation is a second limit and it was put there deliberately. Friction velocity was made to depend on canopy greenness, so the same wind gives a higher friction velocity in July than in February. A threshold estimated over the whole year is therefore a compromise between two different canopies, and the temperature classing that rescues the estimate is doing double duty, separating both the respiration response and the roughness change. On a real site the canopy changes, the fetch changes with wind direction, and instruments are cleaned and replaced. A threshold estimated in one season is not guaranteed to hold in another, which is why Papale et al (2006) estimate the threshold for each quarter of the year and then take the largest of the four. That is a different convention from the single whole-year estimate used here, it will in general return a different number, and neither version is checkable.
The third limit is the one that makes the other two permanent. There is no truth column at a flux tower. The annual figure of -231.46 gC/m2/yr that every method here was scored against exists only because the data were generated. In the field the threshold is chosen by a convention: ninety-nine per cent of the plateau, six temperature classes, twenty bins, the median across classes, the highest of four quarterly values. Change any one of those and the number changes, and there is nothing in the data to say which choice was right. The convention is what makes annual budgets from different towers comparable with each other, and that is a real and sufficient reason to follow it, but comparability is not accuracy.
Where to go next
The record this post leaves behind has 3239 holes in it, they arrive in runs of up to 32 half hours, and they are not missing at random: they are the coldest and quietest part of every night. Filling them is gap filling a flux time series, and the last two sections here are the argument for why that post exists, since the choice of fill changed the annual number by more than the filter did. Once the record is complete, splitting it into gross uptake and respiration is partitioning net flux into GPP and respiration, which fits the same Lloyd and Taylor curve used to generate the data here, to night-time data that has been through this filter. The audit of the finished annual figure, including how much of its uncertainty traces back to the threshold, is checking an annual flux budget.
For the wider machinery: resampling whole nights rather than half hours is the idea behind bootstrapping dependent data: blocks and clusters, and the general treatment of resampled intervals is bootstrap confidence intervals. The non-random hole this filter cuts is a missing-not-at-random problem, whose general form is missing data: MCAR, MAR and MNAR. The litter cluster, doing the same carbon accounting with a balance instead of a sonic anemometer, starts at fitting litter decomposition curves.
References
Goulden ML, Munger JW, Fan SM, Daube BC, Wofsy SC 1996 Global Change Biology 2(3):169-182 (10.1111/j.1365-2486.1996.tb00070.x)
Aubinet M, Grelle A, Ibrom A, Rannik U, Moncrieff J, Foken T, Kowalski AS, Martin PH, Berbigier P, Bernhofer C, Clement R, Elbers J, Granier A, Grunwald T, Morgenstern K, Pilegaard K, Rebmann C, Snijders W, Valentini R, Vesala T 2000 Advances in Ecological Research 30:113-175 (10.1016/S0065-2504(08)60018-5)
Massman WJ, Lee X 2002 Agricultural and Forest Meteorology 113(1-4):121-144 (10.1016/S0168-1923(02)00105-3)
Gu L, Falge EM, Boden T, Baldocchi DD, Black TA, Saleska SR, Suni T, Verma SB, Vesala T, Wofsy SC, Xu L 2005 Agricultural and Forest Meteorology 128(3-4):179-197 (10.1016/j.agrformet.2004.11.006)
Papale D, Reichstein M, Aubinet M, Canfora E, Bernhofer C, Kutsch W, Longdoz B, Rambal S, Valentini R, Vesala T, Yakir D 2006 Biogeosciences 3(4):571-583 (10.5194/bg-3-571-2006)
Richardson AD, Hollinger DY, Burba GG, Davis KJ, Flanagan LB, Katul GG, Munger JW, Ricciuto DM, Stoy PC, Suyker AE, Verma SB, Wofsy SC 2006 Agricultural and Forest Meteorology 136(1-2):1-18 (10.1016/j.agrformet.2006.01.007)
Aubinet M 2008 Ecological Applications 18(6):1368-1378 (10.1890/06-1336.1)
Barr AG, Richardson AD, Hollinger DY, Papale D, Arain MA, Black TA, Bohrer G, Dragoni D, Fischer ML, Gu L, Law BE, Margolis HA, McCaughey JH, Munger JW, Oechel W, Schaeffer K 2013 Agricultural and Forest Meteorology 171-172:31-45 (10.1016/j.agrformet.2012.11.023)
Lloyd J, Taylor JA 1994 Functional Ecology 8(3):315-323 (10.2307/2389824)