While trying to reproduce this validation case with FESTIM I noticed a discrepancy in how the fluxes are applied to D2 (T2). See below for more details.
TLTR: The mass of D2 (T2) is not conserved in this validation case. Which may jeopardise the findings in the published paper.
(I used AI to generate the issue report and figures for the sake of time but I reviewed everything myself)
Bug Description
In test/tests/val-2g/val-2g_trapping.i the dry and wet systems are otherwise
mirror images of each other, but the boundary flux applied to the hydroxyl
variable differs by one term.
flux_on_OT_dry
receives both surface reactions:
expression = '2 * flux_base_on_T2_dry + 2 * flux_base_on_T2O_dry'
flux_on_OT_wet
receives only one:
expression = '2 * flux_base_on_T2O_wet'
flux_on_OT_wet is the only boundary flux on OT_concentration_wet, so the
T2 reaction has no effect on the hydroxyl inventory in the wet system. The same
reaction is nevertheless still active everywhere else in that system:
So in the wet system the reported T2 release is never debited from the solid:
hydrogen leaves the surface without leaving the material. Charge balance breaks
with it, since the electrons are debited but the charged OT⋅ defect is not.
This differs from the model as published
The associated paper (Yang et al., Int. J. Hydrogen Energy 210 (2026)
153551) writes the Q2 surface reaction, Eq. (9), with the hydroxyl balance on
its left-hand side:
½ ∂N(OQ⋅_O)/∂t = −∂N(Q2)/∂t = A ( K₁^Q2 · P_Q2 · C²(Ox_O) − K₋₁^Q2 · C²(OQ⋅_O) · C²(e′) )
i.e. the 2 × rate source on the hydroxyl is part of the equation, exactly as
in Eq. (8) for Q2O, and with no environment-dependent qualifier. The sentence
introducing both equations is explicit that this is not environment-specific:
"Both reactions are reversible on the surface and happen in the dry (with only
hydrogen isotope gas) and wet (with only hydrogen isotope water)
environments."
The only environment dependence stated is the parenthetical one — the
pressures, P_Q2 = 0 in wet and P_Q2O = 0 in dry — which the input already
handles through pressure_T2_wet. That does not remove the hydroxyl sink, which
is driven by the reverse term C²(OQ⋅_O)·C²(e′) and is unaffected by the
pressure being zero.
I could find no assumption in the paper that would justify the omission. The
limitations section lists two model simplifications — a single trapping site,
and a deliberately high initial electron concentration — and neither bears on
this term. The wet T2 flux is also reported and scored against experiment
(40.49% RMSPE in Fig. 4(f)), which is not how an intentionally excluded channel
would be treated.
One observation offered without a claim of causation: the paper's own discussion
of what it cannot explain is about this same channel —
"the release of D2 is faster, with the peak release shifted to lower
temperatures relative to the experiment. This discrepancy cannot be resolved
merely by recalibrating the current material parameters."
— and attributes it tentatively to the electron concentration. Restoring the
missing term does move the wet T2 peak, though it overshoots (to 1173 K against
a measured ~1000 K) rather than resolving the discrepancy, so I would not
present it as the explanation. It may still be worth looking at together.
This can be seen in the repository's own gold file without running anything.
Over the resolved charging window (10–3600 s) val-2g_trapping_calibrated.csv
reports 2.64× more deuterium leaving as T2 than enters as T2O. The dry
system, analysed identically, balances to 5%.
Steps to Reproduce
Runs below used the prebuilt tmap8 2026_07_31 conda package from
https://conda.software.inl.gov/public against the unmodified input at devel
c7643c7.
tmap8-opt -i parameters_trapping_calibrated_validation.params val-2g_trapping.i \
Outputs/file_base=baselinebaseline.csv reproduces gold/val-2g_trapping_calibrated.csv to 1e-11
across every flux and RMSPE column, over an identical 764-step timestep
sequence — i.e. this is exactly the published case.
The patch (conserve_deuterium.diff) makes flux_on_OT_wet an exact mirror
of flux_on_OT_dry. Nothing else is touched:
[flux_on_OT_wet] # OT
type = ADDerivativeParsedMaterial
property_name = 'flux_on_OT_wet'
- material_property_names = 'flux_base_on_T2O_wet'
- expression = '2 * flux_base_on_T2O_wet'
+ material_property_names = 'flux_base_on_T2_wet flux_base_on_T2O_wet'
+ expression = '2 * flux_base_on_T2_wet + 2 * flux_base_on_T2O_wet'
[]Direct measurement of the budget. Adding six pure-diagnostic
postprocessors — the solid inventory ∫(OT + trapped) dx per system, and the
deuterium crossing the surfaces as gas, time-integrated on the solver's own
steps — gives the two sides of the balance. The audited runs are bit-identical
to the un-audited ones (0.00e+00 on every shared column), so the diagnostics do
not perturb the solution. Under conservation the two curves are the same
quantity and must coincide. Worst gap over the run, measured from t = 10 s:
|
dry (control) |
wet |
| as published |
0.05% |
99.05% |
| with the patch |
0.05% |
0.76% |
At the end of charging the wet solid has gained +5.86e10 at/µm² while the gas
balance says deuterium has net left at −9.53e10 at/µm² — opposite signs. Over
the whole run the discrepancy reaches 1.5640e11 at/µm², which matches the
cumulative reported T2 release (both faces, 2 D per molecule) of 1.5636e11 to
0.03%. That pins the mechanism exactly: the missing deuterium is precisely
the T2 that was reported as released.
(The budget is measured from t = 10 s. The initial hydration transient — rate
~5.6e16 over a 1e-4 s step — is not resolved in the output stream, so a
cumulative integral carried through t = 0 is meaningless. The dry control is
what shows the method is sound.)
Impact
The wet agreement in val-2g depends on the missing term. Restoring it:
| RMSPE (final) |
as published |
conserving |
RMSPE_T2_dry |
0.1717 |
0.1717 |
RMSPE_T2O_dry |
0.1774 |
0.1774 |
RMSPE_T2_wet |
0.4171 |
22.694 |
RMSPE_T2O_wet |
0.4597 |
2.2506 |
The dry system is bit-identical between the two runs, as expected — its block
was never touched. In the wet system the T2O desorption peak (3.61e7 at/µm²/s
@ 988 K) disappears entirely; what remains is a steady water-splitting cycle,
T2O absorbed from the 1e-5 Pa residual and T2 emitted, peaking ~185 K above the
measurement.
Physically this is what you would expect: hydration loads far more deuterium
than dry dissociation does, and the T2 recombination rate goes as
C_OT² · C_e², so at high hydroxyl concentration that branch is by far the
cheapest way out of the solid. Once it actually drains the hydroxyl, it drains
it faster than hydration can charge it.
So the patch alone is not a fix for the validation case. It restores
conservation but breaks the agreement, which suggests the calibrated parameter
set — obtained through val-2g_main_PSS_trapping.i, which I have not
inspected — was optimised against the non-conserving wet objective and would
need to be redone. That is the part I can't do from outside, and it is the
reason I am reporting rather than opening a PR.
Two smaller notes for whoever picks this up:
- The conserving variant is stiffer at the end of the ramp, where the T2 sink
drives the hydroxyl toward zero; my runs abort at t = 9395.45 s of 9400 s
(timestep already at or below dtmin). It does not affect the numbers above
— the desorption peak is at 1173 K and the RMSPE window closes at 9165 s —
but the solver settings likely want a look alongside the fix.
- A regression test that the deuterium budget closes in both systems would have
caught this, and would be cheap: the audit postprocessors above are six
blocks and cost nothing measurable.
This appears to be present since the case was added (#357; #387 only touched
the comparison script), and it bears on the val-2g documentation page and on
the associated publication, both of which present dry and wet as one unified
model.
Happy to share the full evidence package — inputs, audited variants, CSVs, and
the figures — or to run any variant you would like checked.
While trying to reproduce this validation case with FESTIM I noticed a discrepancy in how the fluxes are applied to D2 (T2). See below for more details.
TLTR: The mass of D2 (T2) is not conserved in this validation case. Which may jeopardise the findings in the published paper.
(I used AI to generate the issue report and figures for the sake of time but I reviewed everything myself)
Bug Description
In
test/tests/val-2g/val-2g_trapping.ithe dry and wet systems are otherwisemirror images of each other, but the boundary flux applied to the hydroxyl
variable differs by one term.
flux_on_OT_dryreceives both surface reactions:
flux_on_OT_wetreceives only one:
flux_on_OT_wetis the only boundary flux onOT_concentration_wet, so theT2 reaction has no effect on the hydroxyl inventory in the wet system. The same
reaction is nevertheless still active everywhere else in that system:
flux_on_e_wet = 2 * flux_base_on_T2_wetflux_on_OT_wetflux_on_T2_wet = -1 * flux_base_on_T2_wetRMSPE_T2_wetSo in the wet system the reported T2 release is never debited from the solid:
hydrogen leaves the surface without leaving the material. Charge balance breaks
with it, since the electrons are debited but the charged OT⋅ defect is not.
This differs from the model as published
The associated paper (Yang et al., Int. J. Hydrogen Energy 210 (2026)
153551) writes the Q2 surface reaction, Eq. (9), with the hydroxyl balance on
its left-hand side:
i.e. the
2 × ratesource on the hydroxyl is part of the equation, exactly asin Eq. (8) for Q2O, and with no environment-dependent qualifier. The sentence
introducing both equations is explicit that this is not environment-specific:
The only environment dependence stated is the parenthetical one — the
pressures,
P_Q2 = 0in wet andP_Q2O = 0in dry — which the input alreadyhandles through
pressure_T2_wet. That does not remove the hydroxyl sink, whichis driven by the reverse term
C²(OQ⋅_O)·C²(e′)and is unaffected by thepressure being zero.
I could find no assumption in the paper that would justify the omission. The
limitations section lists two model simplifications — a single trapping site,
and a deliberately high initial electron concentration — and neither bears on
this term. The wet T2 flux is also reported and scored against experiment
(40.49% RMSPE in Fig. 4(f)), which is not how an intentionally excluded channel
would be treated.
One observation offered without a claim of causation: the paper's own discussion
of what it cannot explain is about this same channel —
— and attributes it tentatively to the electron concentration. Restoring the
missing term does move the wet T2 peak, though it overshoots (to 1173 K against
a measured ~1000 K) rather than resolving the discrepancy, so I would not
present it as the explanation. It may still be worth looking at together.
This can be seen in the repository's own gold file without running anything.
Over the resolved charging window (10–3600 s)
val-2g_trapping_calibrated.csvreports 2.64× more deuterium leaving as T2 than enters as T2O. The dry
system, analysed identically, balances to 5%.
Steps to Reproduce
Runs below used the prebuilt
tmap8 2026_07_31conda package fromhttps://conda.software.inl.gov/publicagainst the unmodified input at develc7643c7.tmap8-opt -i parameters_trapping_calibrated_validation.params val-2g_trapping.i \ Outputs/file_base=baselinebaseline.csvreproducesgold/val-2g_trapping_calibrated.csvto 1e-11across every flux and RMSPE column, over an identical 764-step timestep
sequence — i.e. this is exactly the published case.
The patch (
conserve_deuterium.diff) makesflux_on_OT_wetan exact mirrorof
flux_on_OT_dry. Nothing else is touched:[flux_on_OT_wet] # OT type = ADDerivativeParsedMaterial property_name = 'flux_on_OT_wet' - material_property_names = 'flux_base_on_T2O_wet' - expression = '2 * flux_base_on_T2O_wet' + material_property_names = 'flux_base_on_T2_wet flux_base_on_T2O_wet' + expression = '2 * flux_base_on_T2_wet + 2 * flux_base_on_T2O_wet' []Direct measurement of the budget. Adding six pure-diagnostic
postprocessors — the solid inventory
∫(OT + trapped) dxper system, and thedeuterium crossing the surfaces as gas, time-integrated on the solver's own
steps — gives the two sides of the balance. The audited runs are bit-identical
to the un-audited ones (0.00e+00 on every shared column), so the diagnostics do
not perturb the solution. Under conservation the two curves are the same
quantity and must coincide. Worst gap over the run, measured from t = 10 s:
At the end of charging the wet solid has gained +5.86e10 at/µm² while the gas
balance says deuterium has net left at −9.53e10 at/µm² — opposite signs. Over
the whole run the discrepancy reaches 1.5640e11 at/µm², which matches the
cumulative reported T2 release (both faces, 2 D per molecule) of 1.5636e11 to
0.03%. That pins the mechanism exactly: the missing deuterium is precisely
the T2 that was reported as released.
(The budget is measured from t = 10 s. The initial hydration transient — rate
~5.6e16 over a 1e-4 s step — is not resolved in the output stream, so a
cumulative integral carried through t = 0 is meaningless. The dry control is
what shows the method is sound.)
Impact
The wet agreement in val-2g depends on the missing term. Restoring it:
RMSPE_T2_dryRMSPE_T2O_dryRMSPE_T2_wetRMSPE_T2O_wetThe dry system is bit-identical between the two runs, as expected — its block
was never touched. In the wet system the T2O desorption peak (3.61e7 at/µm²/s
@ 988 K) disappears entirely; what remains is a steady water-splitting cycle,
T2O absorbed from the 1e-5 Pa residual and T2 emitted, peaking ~185 K above the
measurement.
Physically this is what you would expect: hydration loads far more deuterium
than dry dissociation does, and the T2 recombination rate goes as
C_OT² · C_e², so at high hydroxyl concentration that branch is by far thecheapest way out of the solid. Once it actually drains the hydroxyl, it drains
it faster than hydration can charge it.
So the patch alone is not a fix for the validation case. It restores
conservation but breaks the agreement, which suggests the calibrated parameter
set — obtained through
val-2g_main_PSS_trapping.i, which I have notinspected — was optimised against the non-conserving wet objective and would
need to be redone. That is the part I can't do from outside, and it is the
reason I am reporting rather than opening a PR.
Two smaller notes for whoever picks this up:
drives the hydroxyl toward zero; my runs abort at t = 9395.45 s of 9400 s
(
timestep already at or below dtmin). It does not affect the numbers above— the desorption peak is at 1173 K and the RMSPE window closes at 9165 s —
but the solver settings likely want a look alongside the fix.
caught this, and would be cheap: the audit postprocessors above are six
blocks and cost nothing measurable.
This appears to be present since the case was added (#357; #387 only touched
the comparison script), and it bears on the val-2g documentation page and on
the associated publication, both of which present dry and wet as one unified
model.
Happy to share the full evidence package — inputs, audited variants, CSVs, and
the figures — or to run any variant you would like checked.