Lake Erken — Parity Status and Limitations¶

Status: PARTIAL PARITY (FAIL). The hypsography / inflow / outflow / water-balance is at exact parity. The residual FAIL has two independent causes: (1) a mild, seasonal under-mixing of the summer thermocline — a real but small physical difference affecting the native turbulence/temperature fields — and (2) the selmaprotbas biogeochemistry (BGC) is limited by a FABM library build flag. This page records exactly what is at parity and what is not, with source-level and data-level evidence, so anyone picking up the project knows precisely where it stands.

The case¶

Property	Value
Period / step	1999-02-01 → 2020-12-31 (~22 yr), 3600 s hourly (~192k steps)
Grid	21 m, 42 layers, hypsographic (cross-sectional area varies with depth)
Turbulence	second-order closure (`turb_method=3`, `tke_method=2`, `len_scale_method=8`, `scnd_method=1` quasi-equilibrium, `scnd_coeff=7` = Cheng et al. 2002)
Light	very turbid (`light_extinction` custom: `A=0.58`, `g1=0.5 m`, `g2=1.74 m` — almost all shortwave absorbed in the top ~2 m → a sharp, shallow summer thermocline)
BGC	FABM `selmaprotbas` — 13 interior + 5 benthic state variables
Output	daily means (`time_method: mean`)
FABM numerics	`repair_state: true`, `bottom_everywhere: true`, feedbacks (shade/albedo/drag) off
Lineage	reference built from gotm-lake (older GOTM fork + its own bundled FABM); pyGOTM derives from gotm-model + conda `pyfabm 3.0.0`

Because pyGOTM and the reference come from different lineages, mismatches may originate in gotm-model vs gotm-lake, the translation, the pyFABM coupling, the FABM library build, or the reporting — not every mismatch is a pyGOTM bug.

Variable inventory (current)¶

114 native variables and 41 FABM variables are compared. Native: 60 PASS / 5 MARGINAL / 23 DISCREPANT / 26 BROKEN; FABM: 0 / 0 / 3 / 38. All 64 BROKEN variables are in lake_erken; no other case has any BROKEN.

The 26 native BROKEN variables all trace to one root residual — the summer-thermocline under-mixing described below — either directly or through a sensitive derived/diagnostic quantity. They split as:

Core turbulence/stratification, just over the 0.20 threshold: num (0.245), nuh/nus (0.280), cmue2 (0.238), tke (0.208), as (0.208), SS (0.217), NN/NNT (0.203), G (0.217), avh (0.288), P (0.323).
Sensitive derived diagnostics that amplify the same residual: Rig (0.93 — NN/SS, which blows up wherever SS → 0), the Reynolds-stress diagnostics uu/vv/ww (0.75–0.79), and the SST-driven surface turbulent heat fluxes qe/qh (0.79/0.68) plus the ice fields Hfrazil/Tice_surface.
Momentum: u (0.218), taux/tauy (0.26/0.30), Ekin (0.23).
Mixed-layer-depth diagnostics: mld_bott (0.54), mld_surf (0.45) — see the diagnostic-method note below.

For calibration: temp itself is only DISCREPANT (d_norm = 0.123), rho = 0.141, eps = 0.101. The BROKEN turbulence variables sit just above the 0.20 line because the normalized Fréchet metric is sensitive for these floor-dominated, multi-decade fields; the underlying physical difference is small.

What is at exact parity (verified)¶

Hypsography / water-balance is exact (d_norm = 0.0): Af, Qlayer, Qs, wq, FQ, Qres, zeta, h, salt, and the integrated int_precip / int_evap / int_inflow / int_outflow / int_water_balance terms (the last three PASS within d_norm < 0.003). The momentum equations and bottom friction are hypsography-aware (area-weighted diffusion/advection, wq advection, per-layer bottom drag), mirroring gotm-lake. This is a genuine and significant accomplishment.

The native residual: a mild, seasonal summer-thermocline under-mixing¶

A full-record (8005-day) comparison of the pyGOTM run against the reference NetCDF localises the residual precisely. It is seasonal, not catastrophic:

Season	temp \|err\|	interior `num` ratio (py/ref)	stratification (surface − bottom)
Winter / ice (Nov–Mar)	0.06–0.27 °C	~0.88–1.00	matches
Summer stratified (Jun–Aug)	0.86–1.38 °C	0.44–0.68	py ~2 °C too strong

Annual temp mean-abs-error is 0.55 °C (median 0.35 °C). The largest single point is 6.8 °C, but that is a >99.9th-percentile transient: only ~2 % of points exceed 3 °C and ~0.2 % exceed 5 °C.
The signed bias is a clean dipole: surface +0.49 °C (too warm), bottom −0.20 °C (too cold) on the annual mean; in midsummer the surface runs ~+1 °C and the bottom ~−1 to −1.5 °C. A representative 1-Aug profile: pyGOTM bottom 12.4 °C / surface 22.0 °C (ΔT 9.6 °C) vs reference 13.9 °C / 20.9 °C (ΔT 7.0 °C).
The num deficit is concentrated at the thermocline depth (worst around interface level 32; near-surface and near-bed diffusivity match) and only in the stratified months.
The velocity field matches the reference to ~90 % (near-bottom mean speed: pyGOTM 3.8e-3 vs reference 4.1e-3), and the columns agree in winter and below the thermocline.

Mechanism. During the stratified season pyGOTM holds a sharper, shallower thermocline than the reference: the cross-thermocline turbulent diffusivity is lower, so less heat mixes downward, the surface stays warmer, the bottom stays cooler, and the stratification is ~2 °C too strong. In winter, and below the thermocline, the columns agree. Lake Erken is the case most exposed to this: it is a shallow (21 m), very turbid lake whose shortwave is absorbed almost entirely in the top ~2 m (g2 = 1.74 m), producing an unusually sharp, shallow summer thermocline where any small cross-thermocline mixing difference is maximally amplified. The same closure is only mildly different (MARGINAL/ DISCREPANT) in the deeper, weakly-stratified second-order ocean cases (gotland, nns_annual, resolute).

The second-order closure is identical between the lineages¶

The residual is not a closure-coefficient difference. A direct source diff (gotm-lake/.../turbulence vs gotm-model/code/src/turbulence) shows the closure is the same in both lineages:

Closure component	gotm-lake vs gotm-model
Stability functions `cmue_d.F90`, `cmue_c.F90`	byte-identical
Cheng et al. 2002 coefficients (`CCH02`: `cc1..cc6`, `ct1..ct5`, `ctt`)	byte-identical
`alpha_mnb.F90`, `variances.F90`	identical except Stokes/Langmuir terms (`av`, `aw`, `Px`), which are off for this case
`turbulence.F90` init differences	all in gotm-model-only additions (k-omega namelist, Stokes vectors, renamed enums) that are inactive here

Every closure knob is correctly plumbed from the YAML at runtime (verified through the driver): scnd_method = 1 (quasi-equilibrium), scnd_coeff = 7 (Cheng), length_lim = true, galp = 0.53, compute_c3 = true, Ri_st = 0.25, stab_method = 3. Internal-wave mixing is off in both runs: the case’s iw: block sets only alpha = 0.7 with no method, so iw_model = 0 in pyGOTM, and gotm-lake’s internal_wave.F90 only acts for iw_model == 2. The reference’s interior num shows no ~1e-4 floor, confirming IW mixing is inactive.

Interpretation. With identical closure coefficients and identical, correctly applied configuration, the residual is most consistent with a small numerical difference in how the two codes resolve a sharp, shallow thermocline (advection/diffusion discretisation, accumulated over 22 years) — pyGOTM holds the slightly sharper gradient. There is no evidence-justified code fix: the closure is shared with the 14 passing ocean cases, and forcing pyGOTM to match the reference here would mean adding spurious cross-thermocline mixing, which would degrade those cases. It is recorded as a documented, bounded lineage difference.

Mixed-layer-depth diagnostic: a method-default difference (now configurable)¶

mld_bott and mld_surf are post-hoc diagnostics (diagnostics.F90, mld_method); they do not feed back into the physics. Their d_norm reflects a diagnostic-method default difference between the lineages:

gotm-lake defaults mld_method = 1 (TKE criterion, which also computes the bottom mixed layer); gotm-model (pyGOTM’s basis) defaults mld_method = 2 (critical-Richardson, which computes only the surface MLD and leaves mld_bott = 0). The case YAML overrides neither, so the reference ran method 1 (confirmed: its mld_bott reaches 21 m) while gotm-model’s default is method 2 (mld_bott ≡ 0).
pyGOTM now supports an optional mld: block (method/diff_k/Ri_crit), defaulting to the gotm-model value so the other 22 cases are unchanged. The lake_erken config sets mld: method: 1 to use the same diagnostic as the reference. This is a documented input-config setting (the reference output is untouched); it is diagnostic-only and does not alter any physics field.
Effect (verified): mld_bott d_norm is 0.5423 and mld_surf 0.4477 under method 1 (versus 1.0000 and 0.5950 under the gotm-model default). Both remain BROKEN, because once computed by the reference’s method they faithfully reflect the same summer-thermocline tke residual as the other turbulence variables — they are part of the single coherent story rather than a separate artifact. The case verdict and the 26-native-BROKEN count are unchanged.

FABM coupling fixes already applied (no library rebuild)¶

Two correct, gated fixes live in src/pygotm/fabm/fabm_loop.py (only lake_erken is affected; the other 22 cases are byte-identical):

Daily-mean FABM output (gated on output_reduce_mode): accumulate FABM state and diagnostics over each output window and write the window mean, instead of an instantaneous local-midnight snapshot that biased diurnal quantities.
repair_state clipping (gated on the fabm: repair_state flag): call model.check_state(repair=True) after transport and after the source update, matching gotm-lake do_repair_state. This eliminates a catastrophic selmaprotbas_po phosphate NaN (verified: 0 NaN in the current output).

pyfabm limitations preventing complete BGC validation¶

Two limitations of the stock conda pyfabm library (vs the custom FABM that gotm-lake bundles) mean the selmaprotbas biogeochemistry cannot reach full parity with the reference, independent of the native under-mixing above.

1. variable_bottom_index (benthic coupling). The worst FABM variables are sediment-coupled (selmaprotbas_o2, DO_mg, benthic-regenerated nutrients). gotm-lake’s bottom_everywhere distributes the benthic flux over the whole sloping bed (rhs(k) *= (Af(k)-Af(k-1))/Vc(k) per layer, via set_bottom_index(k)), which requires a FABM library compiled with _FABM_BOTTOM_INDEX_ = -1 (variable bottom index). The conda pyfabm 3.0.0 is built with the default 0; link_bottom_index(...) raises "… compiled without support for variable bottom indices". The reference runs bottom_everywhere only because gotm-lake builds its own FABM with that flag.

2. selmaprotbas alpha/beta phytoplankton parameters (interior BGC). The reference fabm.yaml sets two GOTM-lake-only parameters on each selmaprotbas/phytoplankton instance — alpha (nutrient-uptake half-saturation) and beta (temperature growth-correction). Conda pyfabm does not declare them, so loading the reference YAML verbatim raises an “invalid configuration” error. pyGOTM therefore strips alpha/beta (see src/pygotm/fabm/config.py → _normalized_fabm_config_path) and runs with the upstream defaults, so the interior phytoplankton growth differs from the GOTM-lake-tuned reference. This is why the bundle’s fabm.yaml is a materialized (stripped) copy rather than a byte copy of the reference — the validation runner stages exactly that materialized YAML, and its hash is recorded as fabm_yaml_sha256, so re-running the bundle is reproducible.

Decision: do not rebuild pyfabm locally (conda-forge divergence + maintenance burden on every FABM bump; the variable-bottom build is non-default for the python host and unverified; per-layer evaluation costs nlev × getRates per step; and even with both limitations resolved, full PASS is unlikely given the native under-mixing above). The actionable path is upstream: a conda-forge pyfabm variant that enables variable_bottom_index and exposes the selmaprotbas alpha/beta parameters.

Validation reporting — variable ownership¶

PYGOTM_VARIABLES (the native-vs-FABM section classifier in src/pygotm/validation/tolerances.py) lists the native fields. It mirrors register_all_variables.py plus the runtime lake/air-sea outputs (air-sea, ice, observation, energetics and lake-hypsography fields are native). Normalization is decoupled from ownership — full-range normalization is reserved for the floor-dominated mean-flow/turbulence fields — so section labels and the normalization choice are independent and every reported d_norm is unaffected by the classification.

Reproduction¶

# Re-run lake_erken and regenerate its bundle (~5 min):
conda run -n pygotm python -m pygotm.validation.run_validation --cases lake_erken
# Re-compare existing NetCDFs without re-running:
conda run -n pygotm python -m pygotm.validation.run_validation --all --no-run

Localization diagnostic: compare temp, num, NN at depth between validation/runs/lake_erken/lake_erken.nc and validation/reference/lake_erken/output.nc, grouped by month — the columns track in winter and below the thermocline; the divergence is the summer thermocline.

Bottom line¶

Hypsography / water-balance: exact parity. Native physics: a mild, seasonal under-mixing of the summer thermocline (~0.5 °C annual-mean temp error; stratification ~2 °C too strong in midsummer). The second-order closure source is byte-identical between the lineages and every closure knob is correctly plumbed, so the residual is a small numerical difference in resolving a sharp, shallow thermocline; there is no safe, evidence-justified code fix and the shared closure must not be distorted toward the reference. The mld_bott/mld_surf diagnostics use the reference’s mld_method via the new mld: config block. BGC: two correct FABM-coupling fixes landed; the dominant FABM blocker is bottom_everywhere, which needs a FABM library compiled with variable_bottom_index. The overall case remains FAIL; the physics is close but not at parity.