[Bug]: SPM model SEI on cracks doesn't actually reduce capacity

[lorenzplank]

PyBaMM Version

25.6.1

Python Version

3.12.11

Describe the bug

I am running a coupled degradation simulation based on the https://docs.pybamm.org/en/latest/source/examples/notebooks/models/coupled-degradation.html examples page, but using both DFN and SPM models. My code sets "SEI on cracks": "true" for both models.

Both models report very simular values for "Loss of capacity to negative SEI on cracks [A.h]" and "Loss of capacity to negative SEI [A.h]", so the variable is present and populated for SPM as well as DFN.

However, when I compare the actual capacity loss (time integrated from the charge cycles) to the reported loss contributions:

For DFN, the time integrated capacity loss matches the sum of SEI and SEI on cracks loss.
For SPM, the time integrated capacity loss matches only the SEI loss, and does not seem to include the SEI on cracks loss, even though it is shown as a Loss contribution.

Steps to Reproduce

import pybamm
import matplotlib.pyplot as plt
import numpy as np

# This function integrates the current over time to extract the delivered charge capacity
# for each charge segment (current < 0). Used to compare actual cell capacity loss.
def extract_charge_capacities(solution: pybamm.Solution) -> np.ndarray:
    """
    Extract charge capacities (Ah) from a solution.
    """
    time = solution["Time [s]"].entries
    current = solution["Current [A]"].entries
    discharge_mask = current < 0
    edges = np.where(np.diff(discharge_mask.astype(int)) != 0)[0] + 1
    if discharge_mask[0]:
        edges = np.insert(edges, 0, 0)
    if discharge_mask[-1]:
        edges = np.append(edges, len(current))
    capacities = []
    for i in range(0, len(edges), 2):
        start = edges[i]
        end = edges[i + 1]
        t_seg = time[start:end]
        i_seg = -current[start:end]
        if len(t_seg) > 1:
            capacity_As = np.trapz(i_seg, t_seg)
            capacities.append(capacity_As / 3600)
    return np.asarray(capacities)

# Set up DFN and SPM models with SEI and SEI on cracks enabled.
# This follows the coupled degradation example from the PyBaMM docs.
DFN_model = pybamm.lithium_ion.DFN(
    {
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "particle mechanics": ("swelling and cracking", "swelling only"),
        "SEI on cracks": "true",
    }
)

SPM_model = pybamm.lithium_ion.SPM({
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "particle mechanics": ("swelling and cracking", "swelling only"),
        "SEI on cracks": "true",
    }
)

param = pybamm.ParameterValues("OKane2022")
var_pts = {
    "x_n": 5,  # negative electrode
    "x_s": 5,  # separator
    "x_p": 5,  # positive electrode
    "r_n": 30,  # negative particle
    "r_p": 30,  # positive particle
}

cycle_number = 500
exp = pybamm.Experiment(
    [
        (
            "Discharge at 1C until 2.5 V",  # ageing cycles
            "Charge at 0.3C until 4.2 V",
            "Hold at 4.2 V until C/100",
        )
    ]
    * cycle_number
    + ["Discharge at 0.1C until 2.5 V"],  # final capacity check
)
solver = pybamm.IDAKLUSolver()
DFN_sim = pybamm.Simulation(
    DFN_model, parameter_values=param, experiment=exp, solver=solver, var_pts=var_pts
)
SPM_sim = pybamm.Simulation(
    SPM_model, parameter_values=param, experiment=exp, solver=solver, var_pts=var_pts
)

DFN_sol = DFN_sim.solve()
SPM_sol = SPM_sim.solve()

# Extract loss contributions for both models.
# Note: Both DFN and SPM report nonzero SEI on cracks loss.
DFN_Qt = DFN_sol["Throughput capacity [A.h]"].entries
DFN_Q_SEI = DFN_sol["Loss of capacity to negative SEI [A.h]"].entries
DFN_Q_SEI_cr = DFN_sol["Loss of capacity to negative SEI on cracks [A.h]"].entries
DFN_Q_side = DFN_Q_SEI + DFN_Q_SEI_cr

SPM_Qt = SPM_sol["Throughput capacity [A.h]"].entries
SPM_Q_SEI = SPM_sol["Loss of capacity to negative SEI [A.h]"].entries
SPM_Q_SEI_cr = SPM_sol["Loss of capacity to negative SEI on cracks [A.h]"].entries
SPM_Q_side = SPM_Q_SEI + SPM_Q_SEI_cr

# Plot loss contributions from SEI and SEI on cracks for both models.
plt.figure()
plt.plot(DFN_Qt, DFN_Q_SEI, label="DFN SEI", linestyle="dashed", color="red")
plt.plot(DFN_Qt, DFN_Q_SEI_cr, label="DFN SEI on cracks", linestyle="dashdot", color="blue")
plt.plot(DFN_Qt, DFN_Q_side, label="DFN loss due to side reactions", linestyle="dashdot", color="green")
plt.plot(SPM_Qt, SPM_Q_SEI, label="SPM SEI", linestyle="dotted", color="red")
plt.plot(SPM_Qt, SPM_Q_SEI_cr, label="SPM SEI on cracks", linestyle="dotted", color="blue")
plt.plot(SPM_Qt, SPM_Q_side, label="SPM loss due to side reactions", linestyle="dotted", color="green")
plt.xlabel("Throughput capacity [A.h]")
plt.ylabel("Capacity loss [A.h]")
plt.legend()
plt.show()

# Compare actual integrated charge capacity loss to reported loss contributions.
# For DFN, the actual loss matches the sum of SEI and SEI on cracks.
# For SPM, the actual loss matches only the SEI loss, not the total side reactions.
plt.figure()
plt.plot(-1*(extract_charge_capacities(DFN_sol) - extract_charge_capacities(DFN_sol)[0]), label="DFN Charge capacity loss")
plt.plot(-1*(extract_charge_capacities(SPM_sol) - extract_charge_capacities(SPM_sol)[0]), label="SPM Charge capacity loss")
plt.xlabel("Cycle number")
plt.ylabel("Time integrated Capacity Loss [A.h]")
plt.legend()
plt.show()

Relevant log output