Implementation of SOLPS-ITER code with new Grad–Zhdanov module for D–T mixture

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 nT/nD . 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 nT/nD 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.


Introduction
A deuterium (D) and tritium (T) mixture is used for the fusion reaction in toroidal magnetic devices. Extrinsic seed gases are also used for divertor target heat load control using radiative power dissipation in divertor plasma volume [1,2]. Thus, it is important to be able to simulate D+T+impurity mixtures in fusion machines. It is particularly important to study T distribution in the divertor region, because tungsten (W) sputtering by T is significantly more efficient than by D [3]. Therefore, the T predominance over D near the target plates can lead to the enhanced W sputtering. The difference in the T versus D transport can also affect the T retention pattern with respect to the D retention pattern, which is studied in [4].
The SOLPS-ITER code has been developed for the edge and scrape-off layer transport modelling in axisymmetric toroidal devices, like tokamaks [5,6]. Previously, in the SOLPS-ITER model, a single main light ion and several heavy impurities were assumed [7]. In the present paper, a new Grad-Zhdanov module is implemented into the SOLPS-ITER code. This module is based on the fully multi-ion collisional Grad-Zhdanov closure [8,9]. The kinetic coefficients for parallel ion heat flux, friction and thermal forces and the flow velocity and heat flux dependent parts of the stress-viscosity tensor can be obtained using Grad-Zhdanov theory. Thus, the new Grad-Zhdanov module extends the applicability of the SOLPS-ITER code into arbitrary plasma mixtures. For example, a D+T+Ne mixture can be simulated by the new version of the SOLPS-ITER code.
It should be mentioned that recently a similar ansatz has been applied for test D-T simulations using the SOLEDGE3X code [10]. However, the Grad-Zhdanov stress-viscosity is not presented in that work. Moreover, the multi-temperature generalisation of Grad-Zhdanov method is in progress for the SOLEDGE3X code [11,12].
In the present work, the first D and T test SOLPS-ITER simulations with neon (Ne) impurity in a JET-like configuration are presented, with the aim of showing the capabilities of the new version of the code with the new Grad-Zhdanov module. The D and T spatial separation phenomenon, the contribution from the neutral physics and the Grad-Zhdanov thermal force are studied. Section 2 describes the new Grad-Zhdanov module, with section 3 discussing general simulation parameters and collisional theory applicability. Simulation results are presented in section 4, which is divided into two parts: total plasma transport behaviour (section 4.1) and D and T separation (section 4.2). The contribution from the neutral physics in section 4.2.2 and the contribution from the T ion thermal force with the D ions, together with main ion poloidal flows in section 4.2.3 are separately studied. Finally, in section 5, the possibility for experimental observation of the new effects using H α spectroscopy is considered.

Grad-Zhdanov module
Based on the Grad-Zhdanov multi-species collisional closure [8] and its extension, introduced in [9], the new module has been implemented for the calculation of parallel kinetic coefficients into the 3.0.8 version of the SOLPS-ITER code. In the new module one of the two approaches can be applied [13]: (a) Improved analytical method; (b) Explicit matrix inversion method. When the masses of species are close, as in the present paper for the D+T+Ne mixture, the explicit matrix inversion method should be used. The new Grad-Zhdanov module provides transport coefficients along the magnetic field for: the ion thermal conductivity, transport coefficient for the velocity difference dependent part of the ion heat flux, friction and thermal forces transport coefficients and the coefficients for the flow velocity and heat flux dependent parts of the stress-viscosity tensor [9]. The latter is important to reproduce the neoclassical radial electric field inside the separatix [14]. The new method extends the single-ion implementation, presented in [15], into the multiion case. Moreover, using the new method the single main ion and several heavy impurities assumption was relaxed. Now all equations are written uniformly for all ions, so that the full multi-ion generalisation of the SOLPS-ITER code has been completed. Arbitrary plasma mixtures can now be simulated.

Simulation parameters
The new Grad-Zhdanov module is applied here for the simulation of a D+T+Ne plasma discharge in a JET-like configuration. The magnetic equilibrium and wall geometry are taken from previously published D+Ne SOLPS-ITER simulations for JET-ILW [16], see figure 1. The TRIM-database reflection model for tungsten on the divertor targets and for beryllium on the first wall is applied (figure 1). It should be clear here that the simulation does not intend to match a real JET-ILW experiment, but is rather performed in order to present the new SOLPS-ITER modelling capability for D+T mixtures, particularly in view of the recent DT campaign executed on JET-ILW. Since the SOLPS-ITER JET-ILW simulations reported in [16] were performed by some of the same co-authors featured in the present article, the code set-up and some starting points for the new D+T studies were readily available, making the code runs easier to undertake. Future studies will hopefully be able to focus on detailed comparison between simulations and real experimental results.
The simulation is conducted with drifts (E × B and diamagnetic) and currents enabled [17] and with kinetic neutrals including neutral-neutral collisions [18]. The main plasma parameters are collected in table 1. Power crossing the separatrix is chosen equal to P SOL ≈ 18 MW. Anomalous transport coefficients are taken from [16]. The D and T fluxes across the core boundary are specified, mimicking the source from neutral beam injection (NBI). For these test simulations these fluxes are chosen equal to Γ core D = Γ core T = 7.5 · 10 20 ions s −1 . On JET-ILW two hydrogen isotopes are puffed in different poloidal locations [19]. We intend to exclude the contribution from different puffing locations on the D and T spatial separation. In our test simulation the D and T puffing locations are placed at the same surface above the outer divertor where it is convenient for simulation purposes (cyan zone in figure 1). The puffing fluxes are equal for D and T: Γ puff D = Γ puff T = 1.7 · Figure 1. JET-like modelling geometry including plasma fluid (purple) and kinetic neutral (orange) numerical grids. Beryllium wall-green. Tungsten wall-red. Core boundary-yellow. Puffing surface-cyan. Pumping surfaces-black. Blue surfaces indicate neutral pressure measurements surfaces. 10 22 atoms s −1 . In further studies, it is planned to set puffing locations and puffing flux values according to the exact experimental parameters. The albedo coefficient on the pumping surfaces (black in figure 1) is imposed equal for all species (with the exception of one special case, see section 4.2 and figure 6(d)). As a result, the neutral pressure in the private flux region (PFR) (calculated along the blue surfaces in figure 1) is established at p neut D = 4.7Pa and p neut T = 5.3Pa for D and T respectively. The recycling on all other surfaces is chosen to be 100%. All either D or T species (both ions and atoms), not reflected as atoms from the W divertor surface using the TRIM-database reflection probability, are emitted as part of either D 2 or T 2 molecules, respectively [18,20]. DT molecules are not taken into account in the simulation.
Aiming for a high collisionality regime in the SOL and divertor to keep the collisional theory applicable where the following effects are studied, a relatively high amount of Ne  impurity is seeded in the modelling. The separatrix-averaged Ne concentration is c sep Ne ≈ 2% and Z eff in the confinement region is up to 3 (figure 2). The high value of c sep Ne provides strong divertor radiation, producing low temperature and high density conditions (figure 4), providing the high SOL collisionality.
Following the method in [9] the Knudsen number (Kn) is estimated ( figure 4(a)). As expected, in the confinement region Kn > 1, and the effective mean free-path length of thermal  D ions is larger than the connection length between the top and the bottom of the machine. Note, assuming the poloidal top-bottom connection length is π times larger than the minor radius, where ν * i is defined according to [14]. Thus, Kn > 1/π corresponds to ν * i < ε −3/2 . However, in the SOL, Kn ⩽ 0.1. Therefore, a thermal ion collides 10-100 times on the connection length scale ( figure 4(a)). In the SOL, another criteria can be also used. The Zhdanov parallel heat flux h ∥D + normalized to the convective heat flux moving with the ion thermal velocity √ T i /m D · p D + (figure 4(b)) is the coefficient in the third term in the polynomial correction to the Maxwell distribution function equation (4) in [13]. All coefficients in equation (4) in [13] should be smaller than unity for correct polynomial approximation. In all important parts of the domain the heat-flux-dependent coefficient is found the closest to unity among other coefficients in the distribution function approximation. In the confinement region this criteria is not violated, but the Knudsen number (figure 4(a)) shows that mean-free path is larger than top-bottom connection length, and non-local kinetic effects should be taken into account. Local collisional closure is well applicable in collisional SOL. However, in the regions with high gradient the solution is on the boundary of the collisional theory applicability. When the collisional closure is not applicable neoclassical [21] and Landau damping [22] closures should be correctly applied. Also, one should be careful using such simple analysis. Non-local suprathermal particles can be noncollisional under such conditions and contribute significantly to the heat flux [23], though thermal and friction forces are affected less by non-local suprathermal particles, because these particles weakly interact with the background plasma. The KIPP kinetic code was coupled with 1D version of the SOLPS code for studying contribution from kinetic effects to the transport along the magnetic field [24]. It can be recommended for the important predictive simulation, such as ITER baseline scenario modelling, to choose one flux tube and study if kinetic effects play those major roles.
Note that pronounced detachment occurs on both inner and outer divertor target for this particular choice of simulation parameters. A temperature of the order of 1 eV is observed along the target plates ( figure 3(b)). The peak heat load on the outer target is 2.2 MW m −2 and on the inner target is 1.5 MW m −2 . Qualitatively this simulation corresponds to the Ne 6.0e19 atoms s −1 seeding run in [16] and the 2 MW m −2 peak outer target heat load run in [25]. Note, in contrast to the previous work [16,25] in the present simulation the puffing source is placed above the outer divertor (cyan in figure 1) instead of in the PFR, as it is discussed in this section above. Thus, an additional ionisation source leads to the plasma cooling in the low-field side (LFS) far-SOL ( figure 3(b)).

Stress-viscosity.
Before discussion of the difference in the D and T ion and neutral transport (section 4.2), we consider the impact on the total plasma behaviour of the Grad-Zhdanov closure. Specifically, an importance of the selfconsistent Grad-Zhdanov method application is considered in this section. For this analysis the magnetic flux tube is chosen at the boundary between the near and far SOL (outer midplane (OMP) distance to separatrix is 4 mm; the reference magnetic tube is shown in figure 5(a)), where a significant T i gradient region exists in both inner and outer divertor regions (see figure 5(a)).
Note that in most of the chosen magnetic flux tube the parallel mass-averaged plasma flow velocity (u ∥ ) is around 10% of the sound speed ( figure 5(b)). The velocity is directed towards the inner target from poloidal cell 0 to 64 and towards the outer target from poloidal cell 65 to 95 (figure 5(b)). In the small region near the targets, where a significant recycling ionisation source exists, plasma accelerates up to the sound speed ( figure 5(b)). Also, in the vicinity of the targets recombination exceeds ionization, though most of the particle flux reaches the target. It was found that in the ionisation region the total plasma pressure gradient in the parallel direction is balanced by the ions inertia, while in the SOL further upstream the total plasma pressure is balanced by the parallel component of the ion stress-viscosity tensor, meaning the anisotropic part of the total pressure tensor. Indeed, the multiplication of the Mach and inverse Knudsen numbers (which qualitatively has the meaning of a Reynolds number) estimates what dominates: either the ion inertia or the ion stress-viscosity. In the region where M ≈ 0.1, this multiplication is small in comparison to the region where M ≈ 1 ( figure 5(b)). In the confinement region the radial electric field becomes close to the neoclassical value, and the E×B flux adjusts to zero out the viscousstress tensor. In the SOL, the radial electric field is determined by the sheath electrostatic potential, and the divergence of the ion diamagnetic and E×B flux is balanced by the divergence of the parallel flux, giving rise to non-zero parallel stress-tensor. In this region of the SOL, the parallel pressure gradient is determined by the stress viscosity tensor. The stress-viscosity for each ion D/T/Ne is calculated according to Grad-Zhdanov collisional theory. Thus, a correct evaluation of the viscosity kinetic coefficients influences the plasma density variation on the flux surface, and the drift contribution to the radial transport in the SOL [26]. Therefore, the drift transport across flux surfaces is affected by the Grad-Zhdanov stress-viscosity tensor.
It is important to mention that collisionality is not (figure 5(b)) sufficiently high for the collisional theory applicability everywhere. Note that the stress-viscosity flows are restricted by the sound speed flows of parallel momentum, using so called 'flux limiting coefficients' [27]. In the region of poloidal cell numbers from 46 to 49 and from 55 to 67 the viscosity flux is limited by a factor smaller than 0.9 (with a minimum of 0.6 at poloidal cell 62). In this region the solution is not accurate. Nonetheless, elsewhere the 'flux limiting' does not contribute much and the Grad-Zhdanov viscosity kinetic coefficient is applied.
The poloidal contribution of the Grad-Zhdanov kinematic viscosity b 2 x µ Zh ∥ ≈ 10 3 .. Previously in the SOLPS-ITER code the ion viscosity was calculated based on the modified Braginskii formulation: One can compare the new viscosity coefficient with the one which was calculated using the old method. For example, in the OMP cell 55, for the total plasma viscosity: η  (1). The same result for the low-Z ions is obtained in the chapter 5 in [28]. Note that the transport behaviour discussed in this section is missed in the current version of the SOLEDGE3X code, which accounts only for anomalous cross-field viscosity [10]. Moreover, the parallel stress-viscosity plays a crucial role in the formation of the neoclassical electric field in the confinement region [14].
Besides the strain tensor dependence of the stress tensor, in the confinement region the parallel heat stress is taken into account according to [9]. For comparison with the single ion heat viscosity, one can simply divide : η .00 (taking data at the OMP in the 5th magnetic flux tube inside separatrix as an example). It is close to the simple plasma Pfirsch-Schlüter regime kinetic coefficient, which was previously used in the SOLPS-ITER code [15,17] based on [14]. Comparable Pfirsch-Schlüter regime kinetic coefficients for single ion heat viscosity can be also found in many other works [21,29,30]. Following [15], a similar neoclassical correction α neo is applied in the Grad-Zhdanov module. As previously in the SOLPS-ITER code, the heat stress contribution is only applied on the closed magnetic surfaces. It can be extended to the SOL in the future.
Although the total Grad-Zhdanov stress viscosity for the D+T+Ne mixture is similar to the results for the single main ion (hydrogen isotope) model previously implemented in the SOLPS-ITER code, it is now self-consistently extended for arbitrary plasma mixtures.

Ion heat flux.
According to [9] the total Grad-Zhdanov ion conductivity for the D+T+impurity mixture is supposed to be close to the Braginskii based conductivity, used up until now in the SOLPS-ITER code. In the SOL the difference between the improved ion heat conductivity and the standard SOLPS-ITER ion heat conductivity is observed to be <10%. However, in individual cells the difference can be up to 20%. Inside the separatrix, the difference is ≈17% on the surfaces where Z eff ≈ 3 (figure 2). Thus, the new formulation of the ion heat conductivity does not affect the solution significantly.
In the regions of the large parallel T i gradient, the contribution to the total ion heat flux due to the flow velocity becomes comparable to the conductive ion heat flux. It is the result of the considerable ion flow velocity difference, originating from the balance between the thermal and friction forces. This flux appears due to the D-Ne and T-Ne velocity differences. The contribution to the heat flow due velocity differences between the hydrogen isotopes is small. Overall, the inclusion of D and T as separated species does not result in an significant changes in global plasma behaviour. Some differences in transport with respect to the old SOLPS-ITER are found due to the additional contribution to the ion heat flux and the improved description of thermal and friction forces. The isotope effect contribution to the SOL impurity transport is not a part of this work and left for further studies.

D and T separation
4.2.1. Overview. A predominance of T species over D species is observed in different locations: (a) for ion and neutral species near target surfaces due to smaller T than D effective neutral diffusion and (b) for ion species at the high field side due to the poloidal flows and thermal force effect. It should be mentioned here, the thermal force contribution to the impurity transport is well known [31]. In addition to the total plasma transport (section 4.1) there is different transport behaviour of one hydrogen species with respect to another, which leads to n T /n D ̸ = 1. In figure 6(a) one can see differences between D and T ion densities in SOL and divertor of more than 30%, while their ratio in the confinement region is close to unity. Additionally to the reference 'Realistic ion mass' case (figure 6(a)), other cases in figures 6(b)-(d) are used to support the analysis.
One of the factors that lead to the observed separation of the D and T is a contribution from the thermal force between isotope species and the D/T poloidal flow from the LFS to the HFS. According to [9] as well as to the general nature of the thermal force phenomenon [32], the thermal force is equal to zero for ion species with equal mass. For illustrative purposes (only for the calculation of the kinetic coefficients for thermal and friction forces) we artificially multiply the T mass by 2/3 and make it equal to the D mass ( figure 6(b)). In this case a reduction, with respect to the reference 'realistic ion mass' case ( figure 6(a)), of the relative T to D ion concentration is observed in the SOL, further upstream. Moreover, in the HFS SOL, n T /n D ≈ 1. The opposite limit can be approached by multiplying the T mass by 5 and the Ne mass by 10 for the thermal and friction forces kinetic coefficient calculation ( figure 6(c)). In this case the relative T ion concentration in the SOL, further upstream and in the HFS SOL (figure 6(c)) is even larger with respect to the non-modified case ( figure 6(a)). The Zhdanov-Yushmanov expressions [8], which are based on the heavy impurity assumption (m another ion >> m D ) and used previously in the SOLPS-ITER code, overestimate the thermal force kinetic coefficient and underestimate the friction force kinetic coefficient with respect to the complete Grad-Zhdanov approach with realistic masses [9]. Figure 6(c) demonstrates that the old SOLPS-ITER model with the heavy impurity assumption should not be applied for D-T simulations. The thermal force contribution is discussed in more detail in section 4.2.3.
Note that for all four cases (figure 6) there is a dominance of T ions (and also neutrals, as it is shown in section 4.2.2) within ≈2 cm of both the target plates, with n D > n T being observed further upstream of the plates. This is not, therefore, related to the Grad-Zhdanov thermal force. The effect also appears on much smaller scales than the T i gradient. In addition n T /n D > 1 is found in the PFR for all three cases, in which the same pumping speeds are used (figure 6), whereas n T /n D ≈ 1 in the case of increased effective pumping speed for T species (figure 6(d)). As will be discussed below, these are consequences of the difference in the neutral physics for different hydrogen isotopes.

Neutral physics.
Here we consider the build-up of T ions and neutrals near target surfaces due to the lower effective diffusion of T atoms than for D neutrals. An important first point to make is that the fast atom reflection probability for D and T impinging on the tungsten surface, which is calculated in SOLPS-ITER using the TRIM database [20], is similar for both isotopes. Namely, the average fast atom reflection probability on the divertor plates after ion recycling is ≈0.82 and after atom impact is ≈0.65 (locally the probabilities deviate from these numbers). For technical particle balance concerns in the code, the absorbed ions and atoms are emitted as molecules of the corresponding isotope, i.e. DT molecules are not considered in our simulations. A second simplifying assumption in our model is that the charge exchange (CX) rates are mass-scaled from those for hydrogen, without taking into account effects due to vibrationally-excited states [33].
Considering the neutral particle behaviour following target plate recycling, figure 7 compiles some important neutral parameters along the specific contour directed along flow velocity of D neutral species, i.e. inwards from the inner target as illustrated by the dark blue solid line in figure 6(a). Figure 7(a) gives the estimated ionization and CX frequencies of hydrogen  isotopes along this chosen line, revealing two regions: a diffusion limit for distances <2 cm, where CX dominates over ionization, and a free streaming limit at distances >4 cm where ionization dominates. Another indicator of the two different regimes is the neutral flow velocity (only poloidal and radial components are taken into account) with respect to the sound speed, which is defined here as c neut sa ≡ √ T a /m a ( figure 7(b)). Within the region <2 cm the flow velocity is much smaller than the sound speed, which is an indicator of diffusive flow. Further along the neutral flux line the flow velocity becomes of the order of the sound speed figure 7(b)), which is an indicator of free-streaming flow.
Note that within <2 cm region T dominates over D (ions in figure 6(a) and atoms in figure 7(c)). This is the result of the higher D than T diffusion coefficient. These D and T diffusion coefficients due to the CX process can be estimated as D cx D/T ∼ λ cx D/T c neut sD/T , where λ cx D/T is the D or T mean free-path (i.e. the average length between two CX interactions). Taking into account that λ cx Based on the output from the Eirene kinetic neutral code one can calculate  7(c)). As a result, due to the smaller diffusion coefficient, T atoms are better confined than D atoms in the vicinity of the target plate. Due to the efficient CX process, the condition (4)is also satisfied, and the neutral T predominance over neutral D implies the ion T predominance over ion D. Thus, more T than D particles are involved in the recycling process. Thus, the higher flux of T provides higher ion and neutral densities for T than for D in the case equal D/T flow velocities for neutrals in the region of CX diffusion domination (figures 6(a) and 7(c)).
To illustrate this process a 1D model is proposed, where certain simplifications and assumptions are made. Particle balance, together with the assumption of 100% recycling provide equal neutral and ions fluxes for each isotope: where n a and n neut a are the ion and neutral densities, u ∥ a is the ion parallel velocity (only parallel ion flux is taken into account), the x coordinate is set along the neutral flux line (dark blue solid line in figure 6(a)), the B pol /B represents the projection of the ion parallel velocity on the poloidal plane and the factor α represents the projection of the poloidal ion velocity on the neutral flux line. Note, in the CX diffusion region: due to the large friction between ions (for details see section 4.2.3). Also, the domination of the CX process ensures that: Thus (2) may be written: (5) represents equal diffusion velocities of the D and T species, as can be seen in the simulation in the diffusion region (figures 7(b) and (c)). After integration from the target (target) to the boundary of transition between diffusion and free-streaming regions (boundary): The first '≈' sign in (6) is due to the similar D and T ion fluxes F par T /F par D | boundary ≈ 1 entering the diffusion region (for details see section 4.2.3). Indeed, n neut (3) and (4). For the second '≈' sign in (6) the assumption of constant D neutral flux in the diffusion region is made: Fulfillment of the latter assumption is observed in the simulation. This is a result of the low ionization frequency in this region ( figure 7(a)). Then, the u neut D | target ≈ u diff D and u neut D | boundary ≈ c neut sD | target assumptions are made, where u diff D and c neut sD are respectively the diffusive velocity and neutral sound speed. It is worth to mention that the last '≈' sign in (6) does not appear to be sufficiently accurate: in the simulation c neut sD | target = 1.6u neut D | boundary and u neut D | boundary /u neut D | target = 7 ( figure 7(b)), therefore u diff D ≈ 2.7u neut D | target . However, n neut D | target /n neut D | boundary ≈ 7 is observed in the simulation, which makes the second '≈' sign in (6) reasonable. Finally, n neut T /n neut D | target ≈ 7 0.2 = 1.48 is close to the modelling results ( figure 7(c)). Thus, other assumptions are adequate.
Turning now to distances from the plate >2 cm (figure 7), there is D dominates over T (ions in figure 6(a) and atoms in figure 7(c)). In addition, and as previously noted, there is a transition from the diffusive to the free-streaming regime ( figure 7(b)). During this transition the T atom flux is lower than the D atom flux ( figure 7(c)). Here, the D/T ionization frequencies are close ν ioniz sT /c neut sD < 1. Thus, each T atom has a smaller mean free-path before ionization than a D atom. Therefore, a higher atom D flux is observed >2 cm from target plate. More D than T atoms reach the >2 cm area and ionize there. As a result, n D > n T further from the target in all four figures 6(a)-(d).
To close this subsection, we consider the T ion (and neutral) species prevalence over D ion (and neutral) species in the PFR observed in figures 6(a)-(c). This is a result of the neutral pressure difference between the two isotopes: p neut D = 4.7Pa and p neut T = 5.3Pa (on the blue surfaces in figure 1). The neutral pressure is established according to the throughput and effective pumping speed for each of the species [34]. All complex subdivertor structures can not be included into the 2D SOLPS-ITER model to describe the D and T subdivertor transport up to desired accuracy. Thus, the D and T albedos are (to some extend) free parameters to establish a correct ratio between D and T effective pumping speeds in the PFR, i.e. a correct ratio between D and T PFR neutral pressures for the same D and T throughputs. If we increase the effective T pumping speed by adjusting the T albedo on the pumping surface (black surfaces in figure 1), the T ion (and neutral) species predominance over D ion (and neutral) species in PFR disappears (figure 6(d)). In JET experiment, the p neut D /p neut T ratio can be found using H α spectroscopy [19]. Alternatively, a more complex neutral model including all 3D structures should be applied for D and T subdivertor floating. In the 'increased T pumping' case, n T /n D decreased everywhere in the domain with respect to the 'realistic ion mass' case. Note that n T /n D > 1 is observed upstream and especially at the HFS in the 'increased T pumping' case (figure 6(d)) as also found in the reference case ( figure 6(a)). This effect is the result of the D-T thermal force, which is enhanced in the 'T mass × 5; Ne mass × 10' case (figure 6(c)) and suppressed in the 'T mass × 2/3' case ( figure 6(b)). This phenomenon is discussed in detail in section 4.2.3.

4.2.3.
Thermal force and poloidal flows. As demonstrated, for example, in [25], for the ion ∇B drift directed towards the lower X-point, main ion poloidal flows from the LFS to the HFS (blue arrow in figure 8(a)), are formed (starting from the LFS X-point level-blue circle in figure 8(a)) as a result of drift and recycling flows [25]. Here, we discuss the T build-up in the inner divertor due to both to these poloidal flows and the D-T thermal force. In contrast to the neutral effects discussed in section 4.2.2, the thermal force effect plays its major role on the scale of the parallel T i gradient. The poloidal D and T ion flows in the SOL (blue arrow in figure 8(a)) are fed by the similar D and T fuelling fluxes from the divertor region (red arrow in figure 8(a)), due to the close D and T PFR neutral pressures. Equal D and T puffing and core fluxes also contribute to the D and T poloidal ion fluxes in the SOL. Taking into account that the D poloidal flow stagnation points are found close the T ones in each flux tube, similar D and T poloidal ion fluxes from the LFS to the HFS (with a slight predominance of T flux) are formed. Thus, at the inner divertor entrance the total D and T poloidal ion fluxes averaged from the Xpoint to the outer boundary of the simulation domain relate as: F ion inner T /F ion inner D = 1.15 (poloidal cell 24; figure 8(b)). In the inner divertor the poloidal ion fluxes (magenta arrow in figure 8(a)) increase dramatically fed by the D and T ionisation sources, due to the neutral flux from PFR (green arrow in figure 8(a)). As a result, the ratio of the poloidal ion fluxes for T versus D is F ion T /F ion D ≈ 1.00 − 1.05 in the majority of the inner divertor region: poloidal cell numbers 6-20 in figure 8(b). The similar mechanism equalizes D and T ion fluxes in the outer divertor (poloidal cell numbers 75-92 in figure 8(b)). In the target vicinity the ion fluxes are determined by the local recycling and the T ion fluxes are considerably bigger (details in section 4.2.2).
For the further analysis, the contribution of the perpendicular (drift and anomalous) fluxes should be subtracted from the total ion poloidal fluxes. The remaining parallel contribution into the total ion poloidal flux integrated over the radial coordinate in the divertor is also close for the D and T ( figure 8(b)): F par T /F par D ≈ 1.1 for the most part in the inner divertor. A strong T i parallel gradient is, however, present in the divertor region ( figure 5(a)). This leads to isotope separation due to the thermal force along the magnetic field, i.e. the isotope velocity difference due to the friction force balancing the thermal force between isotopes. We consider this phenomenon in the reference magnetic tube (figure 5(a)). Figure 9(a) plots the thermal force between T and D ions (the thermal forces due to D-Ne and T-Ne interactions are subtracted from the total thermal force), along with the friction force between D and T ions, i.e. the component of the friction force proportional to the D-T parallel flow velocity (u ∥ a ) difference. Except in the target plates vicinity, the thermal force is balanced by the friction force in the inner divertor (similar behaviour can be seen in the outer divertor). It should be noted here that the friction force balancing the thermal force and making ion-ion collisional part of the RHS of the parallel momentum equation close to zero is a known phenomenon observed in the simulation [31]. This thermal force effect leads to the effective impurity leakage from the divertor region. In the area of the high parallel T i gradient ( figure 9(b)) the required ion flow velocity difference between D and T is established to increase the friction force up to the level sufficient to balance the thermal force between the two isotopes. In the inner divertor, in the cold region starting from the target to the 5th poloidal cell, the T i change is small ( figure 9(b)), whereas the collisional frequency is high. As a result of efficient friction force, the ratio u ∥ T /u ∥ D ≈ 1.0 ( figure 9(b)). However, at poloidal cell 6 the T i rises quickly up to the inner divertor entrance (poloidal cell 24). This is associated with a rapid decrease of the u ∥ T with respect to the u ∥ D (both are directed towards the target)(figure 9(b)). As a result, the ion density n T /n D increases, clamping the D and T parallel flow ratio (n T u ∥ T /n D u ∥ D ) close to unity (figure 9(b)), as described previously. This leads directly to n T /n D > 1 on the HFS of the divertor entrance.
This phenomenon can also be described in this way: the total plasma transport with the mass-average flow velocity is defined by the total pressure, the total stress-viscosity tensor, the total inertia and momentum exchange with neutrals. In the frame, which moves with the mass-average flow velocity, D and T defuse with respect to each other. The balance between the thermal and friction forces together with a particle balance defines the D vs T ion density distribution.
The observed T ions predominance over D ions at the separatrix is due to the E × B flux from the HFS where n T /n D > 1( figure 6(a)). It disappears in the 'T mass × 2/3'case ( figure 6(b)).
It is worth mentioning that the T thermal force transport coefficient c (R A T eff) T = 0.28 in case of a 50/50 DT mixture, where the thermal force is written as follows [9]: As shown in [9], in the presence of a non-trace concentration of impurity, the effective T thermal force transport coefficient can be smaller, close to zero or even negative. We also point out that the corresponding Braginskii coefficient for electrons c [32] (in [32] T e is used in (7)). Figure 10 illustrates globally the impact of the Grad-Zhdanov DT thermal force, comparing the distribution of the isotopic ion density ratio (n T /n D ) across the SOL numerical grid of the SOLPS-ITER simulation. n T /n D ≈ 1.3 in the location of the high T i gradient directed against the parallel flow: for example, between the X-point and mid-plane at the HFS ( figure 10(a)). The resulting T build-up at the X-point level on the HFS drives an in-out divertor asymmetry in n T /n D in the 'realistic ion mass' case ( figure 10(a)), whilst when the DT thermal force is suppressed the ratio is more symmetric between inner and outer targets ( figure 10(b)). The smaller n T /n D peak near the OMP (figure 10(a)) is the result of the larger T vs D ion poloidal flux from the LFS to the HFS F ion T /F ion D ≈ 1.17 and a drift driven temperature variations on the flux surfaces. This effect is smaller and we do not study it in details in the present paper.
We close this section with two further important comments related to the ratios of isotopic ion densities. First, the assumption of a large difference in ion masses for the calculation of kinetic coefficients used thus far in the standard SOLPS-ITER model (replicated by the 'T mass × 5, Ne mass × 10' cases in figure 6), leads to an overestimation of the in-out asymmetry of n T /n D (figure 6(c)) -for example, the peak value of n T /n D > 2 at the HFS in these cases. Thus, the old SOLPS-ITER model should not be applied for D-T simulations.
Secondly, n T /n D > 1.5 is observed only in the case in which the D-T thermal force is enabled ('realistic ion mass' vs 'T mass × 2/3'), for p neut T /p neut D ⩾ 1 in the PFR ('Realistic ion mass' vs 'Increased T pumping') and in the region where the neutral physics leads to the T build-up. All three effects contribute comparably to T ions (and neutrals) predominance over D ions (and neutrals) yielding high n T /n D in the vicinity of the inner target plate (figures 6(a) and 10(a)). This could be important for T induced tungsten sputtering and T retention phenomena. It is left for the further studies, when simulations more relevant to the experiment are performed.

Synthetic spectroscopy diagnostic
In JET the isotope ratio averaged along several lines of sight (LoS) passing through the divertor is measured by H α spectroscopy [19]. A similar synthetic diagnostic can be modeled in our JET-like simulations ( figure 11(a)). Here we compute only the total emission of either D or T neutrals without accounting for the specific details of the diagnostic measurement (i.e spectral filter width etc). Nor do the coordinates of the LoS exactly match those in [19]. The intention is only to demonstrate the applicability of such a diagnostic to determine the neutral physics and thermal force impacts on the D-T separation. One could expect similar ratios of D vs T emissions around 656.1 nm lines, when they are considered instead of total radiation, taken in the present work. In our simulations, the major contribution to the collected radiation originates from the highly localized areas where the LoS cross the red-yellow regions ( figure 11(a)). These synthetic measurements are performed for the four simulations discussed in section 4.2.1: 'realistic ion mass' where the unmodified thermal and friction forces are used; 'T mass × 2/3' where the thermal force between D and T is artificially reduced; 'T mass × 5; Ne mass × 10' where the thermal force between D and T is artificially enhanced; 'increased T pumping', where the effective pumping speed for T is increased, and thermal and friction forces are unmodified. Figure 11(b) plots the relative differences between the D and T radiation collected along each LoS. For all cases the T predominance is observed on LoS 1-3 ( figure 11(b)), with the lowest Q D /(Q D + Q T ) observed for the case with the artificially increased D-T thermal force: 'T mass × 5; Ne mass × 10'; the middle Q D /(Q D + Q T ) is found for the cases with the realistic D-T thermal force: 'realistic ion mass' and 'increased T pumping'; the largest Q D /(Q D + Q T ) is seen for the case with the artificially decreased D-T thermal force 'T mass × 2/3'. Thus, the reduction observed on LoS 1-3 is connected to the neutral physics (discussed in subsection 4.2.2), because it is present in all four cases, and enhanced by the D-T thermal force phenomena (discussed in subsection 4.2.3). We note that no significant inner-outer divertor asymmetry is observed in this synthetic diagnostic without the D-T thermal force 'T mass × 2/3'.
The predominance of T over D neutral species along LoS 9-11 is due to the D and T PFR neutral pressures (section 4.2.2). The neutral pressure of each isotope in the PFR is determined by the throughput and effective pumping speed [35]. For the 'increased T pumping' case, there is an increase of Q D /(Q D + Q T ) on all lines with respect to the 'realistic ion mass' case and also in the PFR (LoS 9-11). Thus, for the experimentally determined D-T throughputs, the effective T pumping speed can be adjusted with respect to that of D to reach the correct D-T ratio in the spectroscopy measurements in the PFR, when similar simulations are applied for comparison with real discharges.
We observe finally that at the HFS above the X-point n D /(n D + n T ) ≈ 0.44 for the 'realistic ion mass' case and n D /(n D + n T ) ≈ 0.46 for the 'increased T pumping speed' case. The contribution from the D-T thermal force seems insufficient to be distinguishable by a possible diagnostic measuring ion isotope ratio at the X-point level with respect to the H α spectroscopy measurements.

Conclusions
Using Grad-Zhdanov collisional closure, the complete multiion generalization of the SOLPS-ITER code has been performed and implemented in version 3.0.8 of the code. The assumption of a single main light ion and several heavy impurities is relaxed allowing simulation of a 50/50 D-T mixture in the presence of impurities. The new Grad-Zhdanov module calculates the parallel kinetic coefficients for ion heat flux, friction and thermal forces and the flow velocity and heat flux dependent parts of the stress-viscosity tensor.
Test simulations of the D+T+Ne mixture are conducted for a JET-like configuration, where equal throughputs for D and T species are imposed. The importance of the Grad-Zhdanov stress-viscosity for the total plasma parallel transport is discussed. No significant difference between stress-viscosity kinetic coefficients calculated by the new and old SOLPS-ITER models is found for the particular D+T+Ne mixture. Nevertheless, the stress-viscosity tensor is now calculated selfconsistently for arbitrary plasma mixtures.
Up to 30% T predominance over D is observed (a) in the vicinity of both tungsten targets, (b) in the private flux region, and (c) above the X-point (especially at the HFS), whereas the n T /n D ≈ 1 in the confined region. (a) n T /n D > 1 within <2 cm of both the target plates is observed due to the neutral physics. A simple 1D model is proposed to illustrate the difference in the charge-exchange diffusion for D and T. (b) The dominance of T over D in the private flux region is influenced by the neutral pressure of each hydrogen isotope, which can be equalized by adjusting the T effective pumping speed. (c) The T build-up at the HFS above the X-point is a result of a balance between the D-T thermal and friction forces in the location where the T i gradient in the scrape-off layer is directed against the parallel main ion flow. Furthermore, a higher than 50% T prevalence over D is predicted in the HFS far SOL (near horizontal inner tile) due to the effects listed above.
In the existing SOLPS-ITER model the thermal force was based on the assumption of infinite mass difference between ion species. Approaching the large ion mass ratio limit, a significant overestimation of n T /n D > 2 is found with respect to the complete Grad-Zhdanov closure, where a realistic mass ratio between isotopes is used. The former SOLPS-ITER model is therefore not applicable to D-T simulations.
Finally, a possible identification of the D and T spatial separation by experimental diagnostics is discussed. The contributions of neutral physics and the D-T thermal force effect lead to an n D /(n D + n T ) inner-outer divertor asymmetry, which could be observed by H α spectroscopy. An important direction for future work will be to benchmark these findings against real DT discharges executed in the recent JET DT campaign. At the time of writing this paper, a large number of papers from JET on the results of the DT campaign are in preparation, many of which will provide data suitable to try and model with the new code version.