[Bug]: BasicDFN electrolyte equation does not conserve lithium when the transference number depends on concentration

[repowazdogz-droid]

PyBaMM Version

main at 258fdc8 (2026-09-05); also reproduced on v26.4.1 (d9e5afb) and the PR #5524 head (719763b)

Python Version

3.12.13

Describe the bug

Split out from #5700 at @NicolaCourtier's suggestion.

BasicDFN writes the electrolyte balance as

eps * dc_e/dt = -div(-tor * D_e(c_e,T) * grad(c_e)) + (1 - t_plus(c_e, T)) * a_j / F

(basic_dfn.py, lines 230-233 at 258fdc8), with a purely diffusive flux N_e. That form is equivalent to the conservative cation balance only when t_plus is constant. With t_plus = t_plus(c_e) the conservative form is

eps * dc_e/dt = -div(-tor * D_e * grad(c_e) + t_plus * i_e / F) + a_j / F

and expanding the divergence gives the BasicDFN form plus an extra -i_e * d(t_plus)/dx / F, which BasicDFN omits. Integrating the BasicDFN equation together with the particle equations over the cell gives

dN/dt = -A * integral( t_plus(c_e) * a_j / F ) dx

which is not zero when t_plus varies in x. The modular Full electrolyte submodel keeps t_plus * i_e / F inside the flux (full_diffusion.py, lines 72-76 and 99), which is why the full DFN conserves lithium with the same parameters. #2666 / #2758 fixed exactly this for the modular submodel in 2023; the Basic* models were not part of that change.

Observed (ORegan2022, which has a concentration-dependent t_plus; negative particle diffusivity 3.3e-14; the protocol from #5700; initial SOC 0.45; 20 cells per x-domain, 20 radial shells; IDAKLUSolver(rtol=atol=1e-8)):

Case	max abs(N - N0) / N0
BasicDFN, variable t_plus	2.44e-2
BasicDFN, t_plus fixed at its initial value	1.9e-14
DFN(), same parameters	1.4e-14

• The drift is insensitive to tolerance and mesh: max abs(N - N0) is 6.518e-3 mol at 20/20 for rtol=atol=1e-6, 1e-8 and 1e-10, 6.516e-3 mol at 40/40 and 6.5155e-3 mol at 80/80 (1e-10). Time-integration error is not the explanation, and refinement does not remove it.
• The measured change in N matches the time integral of -A * integral(t_plus * a_j / F) dx from the solution to 1.2e-6 mol out of 6.5e-3 mol (5 s trapezoid).
• Because the constant-t_plus case conserves to round-off with the same node-evaluated source term, this is separate from the finite-volume source-term discretisation question raised in [Bug]: Half-cell and custom models do not conserve lithium #5700. It is a property of the reduced equation, not of how a_j is evaluated on nodes, and a boundary-condition change alone does not remove it (PR Flux boundary conditions in 1D finite volume #5524 leaves the BasicDFN time series byte-identical to main; the model declares no "Flux" boundary condition).

The same reduced form appears in BasicDFN2D (basic_dfn_2d.py line 302), BasicDFNHalfCell (basic_dfn_half_cell.py line 197) and BasicDFNComposite (basic_dfn_composite.py line 323) at 258fdc8. The BasicDFN2D drift @SarahRo reported in #5700 is consistent with this. I have only measured the 1D BasicDFN.

Two reasonable resolutions, and the current state is neither:

1. Put t_plus * i_e / F inside N_e in the Basic* models, as full_diffusion.py does. Electrolyte concentration not conserved if surface form is algebraic #2551 records that div(t_plus * i_e) once hit a concatenation problem in one configuration; I have not tested whether that applies here.
2. State in the Basic* docstrings that they assume a constant t_plus, and warn or raise when a function is supplied. The current docstring cites Marquis 2019 and says nothing about t_plus, and ORegan2022 is accepted silently.

Steps to Reproduce

Run against main 258fdc8; the two lines under "Relevant log output" are what it printed.

import pybamm

experiment = pybamm.Experiment([
    "Rest for 10 minutes (5 seconds period)",
    "Discharge at 1C for 20 minutes (5 seconds period)",
    "Rest for 10 minutes (5 seconds period)",
    "Charge at 1C for 20 minutes (5 seconds period)",
    "Rest for 20 minutes (5 seconds period)",
])

for constant_t_plus in [False, True]:
    model = pybamm.lithium_ion.BasicDFN()
    pv = pybamm.ParameterValues("ORegan2022")
    pv["Negative particle diffusivity [m2.s-1]"] = 3.3e-14
    if constant_t_plus:
        pv["Cation transference number"] = pv["Cation transference number"](
            pv["Initial concentration in electrolyte [mol.m-3]"], pv["Ambient temperature [K]"]
        )
    # lithium per unit area, mol/m2 (BasicDFN has no built-in "Total lithium [mol]")
    n = 0
    for dom, eps in [("Negative", "Negative electrode porosity"),
                     ("Separator", "Separator porosity"),
                     ("Positive", "Positive electrode porosity")]:
        thick = "Separator thickness [m]" if dom == "Separator" else f"{dom} electrode thickness [m]"
        n += pybamm.x_average(model.variables[f"{dom} electrolyte concentration [mol.m-3]"]
                              * pybamm.Parameter(eps)) * pybamm.Parameter(thick)
    for dom in ["Negative", "Positive"]:
        n += pybamm.x_average(pybamm.Parameter(f"{dom} electrode active material volume fraction")
                              * pybamm.r_average(model.variables[f"{dom} particle concentration [mol.m-3]"])) \
             * pybamm.Parameter(f"{dom} electrode thickness [m]")
    model.variables["Lithium per area [mol.m-2]"] = n
    sim = pybamm.Simulation(model, parameter_values=pv, experiment=experiment,
                            solver=pybamm.IDAKLUSolver(rtol=1e-8, atol=1e-8),
                            var_pts={"x_n": 20, "x_s": 20, "x_p": 20, "r_n": 20, "r_p": 20})
    sol = sim.solve(initial_soc=0.45)
    li = sol["Lithium per area [mol.m-2]"].entries
    print(f"constant t+ = {constant_t_plus}: max |N - N0| / N0 = {abs(li - li[0]).max() / li[0]:.3e}")

The refinement table and the comparison with DFN() come from a direct finite-volume sum over cell volumes rather than x_average; code and per-run outputs are in the package linked from #5700 (https://github.com/repowazdogz-droid/pybamm-5700-lithium-inventory, commit adaa877).

Relevant log output

constant t+ = False: max |N - N0| / N0 = 2.442e-02
constant t+ = True: max |N - N0| / N0 = 1.982e-14

[rtimms]

Thanks, I would prefer option 1.