McNabb–Foster Equations in Hydrogen Diffusion

• McNabb–Foster equations are a generalized set of kinetic equations modeling hydrogen diffusion with reversible saturable traps in metals.
• They employ both finite-element and finite-difference methods to simulate coupled lattice diffusion, trap occupancy, and experimental permeation and desorption processes.
• Applications in ferritic steels highlight the critical interplay of trap density and binding energies, stressing the need for uncertainty quantification in inverse parameter estimation.

The McNabb–Foster equations are a generalized set of equations for hydrogen diffusion in the presence of saturable traps. In contemporary arXiv treatments of hydrogen transport in steels, they provide the governing framework for coupled lattice diffusion, trap occupancy evolution, permeation, room-temperature desorption, and thermal desorption spectroscopy (TDS). The framework is explicitly kinetic rather than purely local-equilibrium, and recent studies emphasize both its forward-modeling power and a central inverse-problem limitation: electrochemical permeation (EP) and TDS data can often be reproduced by multiple combinations of trap density and trapping parameters, so individual parameters are not necessarily uniquely identifiable (Winzer, 21 Jul 2025, Nunzio, 24 Apr 2025).

1. Governing equations and conservation structure

A recent finite-element treatment describes the McNabb–Foster framework as the “generalised set of equations for hydrogen diffusion in the presence of saturable traps” (Winzer, 21 Jul 2025). In that formulation, the lattice hydrogen concentration $C_L$ evolves according to

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$

while each trap family $i$ satisfies

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$

Here $\theta_i=C_i/N_i$ is the fractional occupancy of trap family $i$, $N_i$ is the trap-site density, $K_i$ is the trapping rate constant, and $\lambda_i$ is the detrapping rate constant. The total hydrogen concentration may be written as

$C_{\text{tot}} = C_L + \sum_{i=1}^{n} N_i\theta_i,$

so the lattice equation is equivalently a conservation statement,

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$0

In fitting examples for ferritic steels containing $$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$1 particles, the general system was reduced to a single reversible saturable trap,

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$2

This structure distinguishes the McNabb–Foster description from a purely Fickian effective-diffusivity treatment. It retains explicit exchange between mobile lattice hydrogen and trapped hydrogen, and therefore represents transient storage and release rather than absorbing all trap effects into a single apparent diffusivity (Winzer, 21 Jul 2025).

2. Kinetic parameters, occupancies, and related formulations

The same recent treatment gives Arrhenius forms for the trapping and detrapping rate constants,

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$3

with trap binding energy

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$4

In later notation used in the fitting discussion, this becomes

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$5

The lattice diffusivity was fixed as

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$6

with

$$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$7

These constants were used in the finite-element model for ferritic steels (Winzer, 21 Jul 2025).

A related arXiv formulation presents the same physical content in terms of occupancies of lattice and trap sites rather than directly in terms of $$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$8 and $$\frac{\partial C_L}{\partial t} = D_L \nabla^2 C_L - \sum_{i=1}^{n} N_i \frac{\partial \theta_i}{\partial t},$$9 (Nunzio, 24 Apr 2025). There,

$i$0

and

$i$1

The transport equation is written as

$i$2

In that specialization, adsorption and desorption are occupancy-dependent and thermally activated. The trapping term carries the barrier $i$3, whereas the detrapping term carries the barrier $i$4. The same paper states that diffusion and adsorption have the same activation barrier $i$5, while desorption must overcome both diffusion and trap binding. This asymmetry is central to its representation of slow release from deep traps during TDS (Nunzio, 24 Apr 2025).

The relation to Oriani’s local-equilibrium approximation is explicit. One paper states that Oriani’s simplification retains only $i$6, so the separate barriers $i$7 and $i$8 cannot be extracted, and therefore deliberately uses the full McNabb–Foster kinetic model to “capture as much detail as possible regarding the trap energetics” (Winzer, 21 Jul 2025). Another presents the Oriani occupancy relation

$i$9

but keeps finite trapping and detrapping kinetics in its numerical model (Nunzio, 24 Apr 2025). This suggests that the McNabb–Foster framework is treated as the more general kinetic description, with Oriani recovered only in appropriate fast-equilibrium limits.

3. Boundary-value problems for permeation, desorption, and TDS

The McNabb–Foster equations are used as a unified framework for EP and TDS. In the finite-element implementation applied to ferritic steels, the EP problem used the initial condition

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$0

a fixed charging-side lattice concentration at the cathodic surface,

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$1

and

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$2

at the anodic side. The observable current density at the exit surface was computed from the flux $$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$3 as

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$4

The same study notes that a simple Fickian baseline,

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$5

is inadequate for the measured curves, which deviate significantly from the analytical Fickian form (Winzer, 21 Jul 2025).

For TDS, the same paper imposed zero lattice concentration on both surfaces,

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$6

took the initial total hydrogen concentration from the integrated TDS spectrum, and assumed that all initial hydrogen is trapped because any lattice hydrogen would have escaped before heating started. The measured temperature history was imposed experimentally, and the total desorption signal was taken as the sum of the fluxes from both surfaces (Winzer, 21 Jul 2025).

A related sensitivity-analysis study uses the same unified sequence of charging/permeation, free desorption, and TDS. During permeation it imposes

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$7

and during free desorption and thermal desorption both surfaces are treated as zero-concentration sinks. The surface flux is computed as

$$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$8

That paper emphasizes that a TDS peak is not a pure trap-emptying event but arises from coupled detrapping, lattice diffusion, specimen thickness, and heating rate (Nunzio, 24 Apr 2025).

4. Numerical realization and parameter-estimation workflows

One arXiv implementation solves the McNabb–Foster equations in FESTIM, using a 1D through-thickness membrane geometry discretized with 100 unidimensional linear elements of equal length. Typical time steps were 100 s for EP and 5 s for TDS. In all fitting examples, a single type of trap was assumed. Least-squares fitting was performed using lmfit, based on the Levenberg–Marquardt algorithm, and the optimized parameters were usually $$\frac{\partial \theta_i}{\partial t} = K_i C_L (1-\theta_i) - \lambda_i \theta_i.$$9, $\theta_i=C_i/N_i$0, and $\theta_i=C_i/N_i$1 (Winzer, 21 Jul 2025).

That same study also documents what was not explicitly specified: no objective-function equation was printed, no weighting strategy was described, no parameter bounds were reported, no initial guesses were reported, and no convergence criteria were reported. Instead, the practical strategy for identifiability assessment was to fix one parameter—especially $\theta_i=C_i/N_i$2—over broad ranges and re-optimize the others, then examine whether good fits persisted (Winzer, 21 Jul 2025).

A second implementation uses an explicit finite-difference scheme rather than finite elements. The specimen thickness is divided into $\theta_i=C_i/N_i$3 uniformly spaced nodes,

$\theta_i=C_i/N_i$4

and the inner-node update is

$\theta_i=C_i/N_i$5

with

$\theta_i=C_i/N_i$6

The time step must satisfy the von Neumann condition

$\theta_i=C_i/N_i$7

and additional local mass-balance constraints are enforced so that trapping and detrapping increments do not exceed the available lattice or trap populations (Nunzio, 24 Apr 2025).

Taken together, these implementations show that the McNabb–Foster equations are not tied to a single numerical method. What remains invariant is the coupled diffusion–trapping–detrapping structure and the use of the same governing system across permeation and desorption protocols.

5. Identifiability, parameter coupling, and non-uniqueness

The most important recent result concerning the McNabb–Foster equations is not the existence of a forward solver but the non-uniqueness of inverse parameter identification. In ferritic steels containing $\theta_i=C_i/N_i$8 particles, EP and TDS curves were fitted with multiple combinations of trap density and energetic parameters, leading to the conclusion that the system was overdetermined and that it was not possible to determine the individual trapping parameters using this procedure (Winzer, 21 Jul 2025).

The EP example for specimen T3 is explicit. With

$\theta_i=C_i/N_i$9

the fitted values were

$i$0

With

$i$1

the fitted values were

$i$2

Across all EP measurements,

$i$3

For TDS, good fits were obtained for specimen T4 over the much wider range

$i$4

Within this range, $i$5 stayed relatively consistent, $i$6 increased approximately linearly with $i$7, and $i$8 decreased linearly with $i$9. The gradient of $N_i$0 versus $N_i$1 was $N_i$2 eV for the T4 example and $N_i$3 to $N_i$4 eV across materials (Winzer, 21 Jul 2025).

The mathematical explanation is the coupling

$N_i$5

which, after substituting the Arrhenius forms, becomes

$N_i$6

At $N_i$7 K, $N_i$8 eV; at TDS peak temperatures $N_i$9–$K_i$0 K, $K_i$1 to $K_i$2 eV. These values are consistent with the observed slopes. A plausible implication is that EP and TDS constrain a combined trapping response more strongly than they constrain the individual quantities $K_i$3, $K_i$4, and $K_i$5 (Winzer, 21 Jul 2025).

The comparison with Kissinger analysis reinforces this point. One study applied

$K_i$6

and, after peak deconvolution, reported values roughly in the range $K_i$7 to $K_i$8 eV, with most values around $K_i$9 to $\lambda_i$0 eV; its abstract summarizes this as a trap binding energy of about $\lambda_i$1 eV, albeit with a high degree of uncertainty (Winzer, 21 Jul 2025). A separate sensitivity analysis also states that TDS peak temperature depends on trap binding energy, trap density, specimen thickness, and heating rate, so peak position alone is not a unique measure of binding energy (Nunzio, 24 Apr 2025).

6. Scope, applications, and terminological boundaries

Within the recent arXiv literature, the McNabb–Foster equations are applied to hydrogen transport in ferritic steels and are used specifically to interpret EP and TDS in systems containing microstructural traps. In ferritic steels containing $\lambda_i$2 particles, measurements showed that fine particles < 5 nm slowed hydrogen diffusion significantly, whereas coarser particles with average diameter > 10 nm had little or no effect. The McNabb–Foster framework reproduced these observations in a forward sense but did not isolate a unique trap parameter set in the inverse sense (Winzer, 21 Jul 2025).

The framework also carries a clear methodological warning. One study explicitly recommends against over-interpreting a single “best-fit” McNabb–Foster parameter set from EP/TDS alone and instead suggests reporting admissible parameter ranges, correlations between parameters, and uncertainty or non-uniqueness (Winzer, 21 Jul 2025). Another reaches a related conclusion from sensitivity analysis: high-energy trap binding energy can strongly affect TDS while having little effect on permeation time lag once traps are saturated, and TDS peak temperature is influenced by geometry and trap density as well as binding energy (Nunzio, 24 Apr 2025).

A common terminological confusion arises from the word “Foster.” The phrase “McNabb–Foster equations” does not appear in the paper on the Foster–Hart measure of riskiness, where the relevant equation is

$\lambda_i$3

and it also does not appear in papers on non-Foster electromagnetic media, temporal metastructures, or photonic time crystals (Riedel et al., 2013, Pacheco-Peña et al., 2023, Li et al., 31 Aug 2025). In materials science usage, by contrast, the term refers to the kinetic trapping–diffusion framework for hydrogen in metals.

In that restricted and technically specific sense, the McNabb–Foster equations denote a non-equilibrium, reversible, saturable trap model for hydrogen transport. Their principal value lies in unifying diffusion, trapping, detrapping, permeation, and desorption within one constitutive system. Their principal limitation, as recent arXiv studies emphasize, is that good agreement with EP and TDS does not by itself establish unique microscopic trap parameters.