postprocessing.py:182:
concentration_mean = (ds["concentration"] * ds["thickness"]).sum(dim="layer") / ds["thickness"].sum(dim="layer")
xarray's .sum() skips NaN by default, so a layer whose concentration is NaN (inactive/dry cell) contributes nothing to the numerator, but its thickness still counts in full toward the denominator. The thickness-weighted mean is therefore diluted toward 0 wherever any layer is NaN, and the size of the error scales with that layer's thickness.
Compare nlmod.layers.aggregate_by_weighted_mean_to_ds, which masks the denominator with s_thk.where(~np.isnan(source_ds[var_name])) — the same weighting done correctly.
Likely fix: mask the denominator by ds["concentration"].notnull(), or use conc_filled deliberately if inheriting the value from below is intended.
Per the repo review doctrine, a regression test reproducing the bias against the current baseline should be written before the fix lands. Found while writing the nhflotools test suite; deliberately not pinned by a test so the current behaviour is not entrenched.
postprocessing.py:182:xarray's
.sum()skips NaN by default, so a layer whose concentration is NaN (inactive/dry cell) contributes nothing to the numerator, but its thickness still counts in full toward the denominator. The thickness-weighted mean is therefore diluted toward 0 wherever any layer is NaN, and the size of the error scales with that layer's thickness.Compare
nlmod.layers.aggregate_by_weighted_mean_to_ds, which masks the denominator withs_thk.where(~np.isnan(source_ds[var_name]))— the same weighting done correctly.Likely fix: mask the denominator by
ds["concentration"].notnull(), or useconc_filleddeliberately if inheriting the value from below is intended.Per the repo review doctrine, a regression test reproducing the bias against the current baseline should be written before the fix lands. Found while writing the nhflotools test suite; deliberately not pinned by a test so the current behaviour is not entrenched.