A new Grad–Zhdanov module is implemented into the SOLPS-ITER code for calculation of the parallel kinetic coefficients. A complete, multi-ion generalization is performed, relaxing the heavy impurity assumption. A JET-like D+T+Ne test simulation is conducted to demonstrate the ability to model a 50/50 deuterium (D) and tritium (T) mixture in SOLPS-ITER. More than 30% T build-up with respect to D is observed in different parts of the simulation domain, in particular at the high field side. A T predominance over D is also observed near both the inner and outer targets. It is a result of the different effective diffusion, dominated by charge-exchange processes, for the D and T neutral species. A simple 1D model is proposed to describe this phenomenon. The contribution to the differing D and T distribution in different regions from the Grad–Zhdanov thermal force is also studied. Due to the thermal force and the D/T poloidal flow from the low field side to the high field side, the prevalence of the T species over D is found at the high field side at the X-point level. The latter leads to an inner–outer divertor asymmetry of . However, the effect is relatively small due to the close D and T masses. It is further shown that the infinite ion mass difference limit, which is used for the derivation of the Zhdanov-Yushmanov analytical expressions and applied for the thermal force coefficient calculation previously used in the SOLPS-ITER code, overestimates significantly. Thus, the old SOLPS-ITER model should not be applied for D–T simulations. Finally, possible experimental studies of the D and T spatial separation due to the new effects revealed by this modelling are discussed.