CoverageDependentSurfPhase does not satisfy Gibbs-Duhem relation

[speth]

Problem description

The coverage-dependent-surface thermo model adds coverage-dependent enthalpy and entropy corrections to each surface species independently. The resulting chemical potentials do not satisfy the Gibbs-Duhem relation at constant T and P,

$$\sum_k X_k d\mu_k = 0,$$

because the per-species corrections lack the cross-species coupling required for the $\mu_k$ to be partial molar Gibbs energies of a single $G(T, P, \mathbf{N})$. The violation is large — the residual is the same order of magnitude as the chemical-potential changes themselves (i.e. essentially no cancellation). Conceptually, I think this is similar to the issue in #2132 where per-species rather than binary interaction coefficients introduce problems.

Steps to reproduce

import cantera as ct

# Same species and same file; "ideal" is a plain ideal-surface (consistent
# control), "covdep" adds the coverage-dependent corrections.
for phase_name in ["ideal", "covdep"]:
    s = ct.Interface("test/data/copt_covdepsurf_example.yaml", phase_name)
    s.TP = 700.0, ct.one_atm
    s.coverages = {"OC_Pt": 0.1, "CO2_Pt": 0.5, "C_Pt": 0.1, "O_Pt": 0.3}
    X0 = s.X
    print(f"\n=== {phase_name} ===")
    active = [k for k in range(s.n_species) if X0[k] > 1e-6]
    i = active[0]
    for j in active[1:]:
        dx = 1e-4 * min(X0[i], X0[j])
        Xp = X0.copy(); Xp[i] -= dx; Xp[j] += dx
        Xm = X0.copy(); Xm[i] += dx; Xm[j] -= dx
        s.TP = 700.0, ct.one_atm; s.coverages = Xp; mu_p = s.chemical_potentials
        s.TP = 700.0, ct.one_atm; s.coverages = Xm; mu_m = s.chemical_potentials
        gd = sum(X0[k]*(mu_p[k]-mu_m[k]) for k in range(s.n_species))
        scale = max(abs(X0[k]*(mu_p[k]-mu_m[k])) for k in range(s.n_species))
        print(f"  {s.species_names[i]} <-> {s.species_names[j]:7s}: "
              f"sum_k X_k dmu_k = {gd: .4e} J/kmol   (rel {gd/scale: .2e})")

Behavior

=== ideal ===
  OC_Pt <-> CO2_Pt : sum_k X_k dmu_k = -3.4571e-07 J/kmol   (rel -2.97e-09)
  OC_Pt <-> C_Pt   : sum_k X_k dmu_k = -1.1176e-09 J/kmol   (rel -9.60e-12)
  OC_Pt <-> O_Pt   : sum_k X_k dmu_k = -3.4571e-07 J/kmol   (rel -2.97e-09)

=== covdep ===
  OC_Pt <-> CO2_Pt : sum_k X_k dmu_k =  2.8764e+05 J/kmol   (rel  1.09e+00)
  OC_Pt <-> C_Pt   : sum_k X_k dmu_k = -1.0002e+05 J/kmol   (rel -1.64e+00)
  OC_Pt <-> O_Pt   : sum_k X_k dmu_k = -7.0237e+04 J/kmol   (rel -2.07e+00)

The ideal-surface control satisfies Gibbs-Duhem to round-off; the coverage-dependent model's relative residual is order 1. Expected: a residual approaching zero.

System information

• Cantera version: 4.0.0a1 / main branch at

Additional context

Found by a new gibbs_duhem_const_T_P consistency test; currently marked known-failure there. A consistent formulation would require the coverage-dependent corrections to be derived from a single coverage-dependent excess Gibbs energy.

Issue draft written by Claude Opus 4.8 after extensive discussion.