Coupled KIPP-EDGE2D modelling of parallel transport in the SOL and divertor of inter-ELM JET high radiative H-mode plasma

The Kinetic Code for Plasma Periphery (KIPP) models parallel (along magnetic field lines) propagation of charged particles in the scrape-off layer (SOL) and divertor of tokamaks. An iterative coupling between KIPP and a 2D edge fluid code EDGE2D, which in turn is coupled to the Monte-Carlo solver EIRENE for neutrals, was used to achieve a converged KIPP-EDGE2D-EIRENE solution. The original EDGE2D-EIRENE solution simulated SOL and divertor of JET high radiative inter-edge localized mode H-mode plasma conditions with strong nitrogen injection, leading to partial detachment at divertor targets. This work is a continuation of earlier studies of modelling kinetic electrons (Chankin et al 2018 Plasma Phys. Control. Fusion 60 115011) and ions (Chankin et al 2020 Plasma Phys. Control. Fusion 62 105022) with KIPP. For numerical reasons caused by large cell-to-cell plasma parameter variations near entrances to divertors, multipliers for parallel electron and ion conductive power fluxes (KIPP/EDGE2D ratios) which are passed onto EDGE2D, could only be used in the main SOL, outside divertors. There, the heat flux limiting effect led to an increase in maximum plasma temperatures in the main SOL and a decrease in power fluxes to divertor targets. Results of the coupling studies are consistent with earlier studies, suggesting that under investigated JET plasma conditions kinetic effects of charged particle parallel propagation do not drastically change target power deposition at divertor targets calculated by EDGE2D-EIRENE along.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The Kinetic Code for Plasma Periphery (KIPP) is a continuum Vlasov-Fokker-Planck code for parallel plasma transport (along magnetic field lines) in the scrape-off layer (SOL) and divertor. The main equations, including the normalization of parameters, are described in [1]. The code combines an implicit 2nd order scheme for a full non-linear Coulomb collision operator with an explicit 2nd order scheme for the parallel free-streaming. Results of the code benchmarking can be found in [2] and references therein. In the present paper the code was run in the 1D2V mode, with one spatial coordinate (along magnetic field lines) and two velocity variables: parallel and gyro-averaged perpendicular velocities.
The present paper is a continuation of two earlier studies on the modelling kinetic electrons [3] and ions [4] with KIPP. In both studies the EDGE2D-EIRENE solution simulating SOL and divertor of JET high radiative H-mode inter-edge localized mode (ELM) plasma conditions with 8 MW of input power into the computational grid and strong nitrogen (N 2 ) injection, leading to the partial detachment at divertor targets [5], was used as a starting point for KIPP runs. Distribution functions ( f e and f i ), analysed in [3,4], reveal that non-Maxwellian features range from high energy tails in the main SOL to bump-on-tail features owing to strongly non-local transport, super-imposed onto 'Maxwellian cores' in the divertor under detached conditions. The latter are caused by very high Maxwellization rates for the bulk of charged particles under high density and low temperature conditions. Here only macroscopic plasma parameters obtained from EDGE2D-EIRENE solutions and extracted from KIPP runs are presented. By 'plasma temperature' in KIPP runs one understands 2/3's of energy content ⟨ f e,i v 2 e,i /2⟩ (with ⟨. . .⟩ denoting averaging over the distribution function) divided by density n e,i in the reference frame moving with the plasma species at velocity ⟨v ||e,i ⟩.
Here the aim was to apply an iterative coupling scheme successfully tested earlier (KIPP-scrape-off layer plasma simulation (SOLPS) coupling, see [6]) to JET radiative divertor conditions. In [6] only ∼5 iterative steps, with KIPP supplying SOLPS with kinetic transport coefficients and SOLPS returning to KIPP plasma parameter profiles maintained by sources, were sufficient to reach a steady state coupled KIPP-SOLPS solution. In the present work, due to numerical problems caused ultimately by very large cell-to-cell parameter variations near entrances to, and inside divertors (see details in section 3), applying such an iterative coupling scheme in divertor regions was impossible. At the same time, the scheme worked well in the main SOL, outside of divertors, where conductive electron and ion power fluxes are known to be lower than those calculated by fluid codes ('heat flux limiting').
In this paper, the setup of EDGE2D-EIRENE cases and the KIPP-EDGE2D-EIRENE (referred to as KIPP_EDGE2D) iterative coupling algorithm is described in section 2. Challenges to the implementation of the iterative KIPP-EDGE2D coupling scheme in divertor regions are described in section 3. Parallel power flux and temperature profiles are presented in section 4. Power balance in EDGE2D solutions corresponding to progressing KIPP-EDGE2D iterations is considered in section 5. The results of the work are summarized in section 6. Figure 1 shows an expanded view of the EDGE2D grid in the divertor region, indicating numbered cells: only odd numbers out of total 16 cells in the radial direction are indicated. These cells belong to 'poloidal rings' in the SOL, labelled s01 to s16, according to EDGE2D nomenclature. Arrows indicate the direction of counting parallel distances, from the inner to outer target. Note that such a counting is opposite to that adopted in EDGE2D (from the outer to inner target). The KIPP grid is based on the EDGE2D grid. However, for the sake of obtaining smoother profiles in KIPP solutions, each EDGE2D cell is divided into two in the poloidal/parallel direction, with the corresponding interpolation for plasma parameters onto the new, finer grid. The number of cells in KIPP is therefore doubles that of EDGE2D. The reader is referred to [3,4] for details of KIPP kinetic calculations along magnetic field lines.

Setup of EDGE2D-EIRENE cases and KIPP-EDGE2D iterative coupling
The latest EDGE2D version allows the introduction of 2D arrays specifying KIPP/EDGE2D multipliers for parallel transport coefficients for electrons and ions, different in each cell. As pointed out above, in this work the multipliers were applied only to EDGE2D transport coefficients for cells in the main SOL. The transport coefficients affected by the multipliers were heat conductivities for electron and main ion (deuterium) parallel conductive power fluxes. As explained in [4], impurities (seven nitrogen and four beryllium species) were not treated kinetically, their distribution functions were assumed to be drifting Maxwellian. They were accounted for in the calculation of Rosenbluth potentials, which in turn were used to calculate transport coefficients is velocity space for electrons and main ions. The multipliers were calculated by dividing conductive power fluxes obtained in KIPP runs by those obtained in preceding EDGE2D-EIRENE runs.
Near stagnation points in the main SOL, where electron and ion temperatures, T e and T i , reached maxima, the multipliers could have negative numbers for some cells. This can be explained by differences in calculations of power fluxes in EDGE2D and KIPP. In EDGE2D analytical formulas for electron and main ion conductive power flux coefficients χ e and χ i , where used, yielding conductive power fluxes in the direction opposite to parallel T e and T i gradients. In KIPP due to non-local effects of the power transport, T e,i maxima did not necessarily coincide with locations of zero conductive power fluxes. Therefore, in some cells the multipliers had negative values. It was, however, found in EDGE2D runs that negative χ e and χ i values quickly lead to instabilities and failure of cases. For this reason negative multipliers were not applied to such cells, that is, they were set to 1. Away from stagnation points, where strong parallel T e and T i gradients develop, the calculated multipliers were always positive and therefore were fully applied. The consequence of effectively not applying the multipliers to some cells near stagnation points can be seen  in figures comparing EDGE2D and KIPP conductive power profiles in section 4, showing a significant difference between these profiles.
The first EDGE2D case, with all multipliers set to 1, served as the 0th iteration EDGE2D case, which is essentially the original case in [5], but continued under the newer EDGE2D and EIRENE versions. Based on this, the KIPP run generates KIPP/EDGE2D multipliers for the next, 1st iteration EDGE2D run, and so on. As in [6], only five KIPP-EDGE2D iterations were sufficient to reach converged solutions. Figure 2 shows radial profiles of T e , T i and electron density n e across cells at the outer midplane (OMP) for the 0th iteration EDGE2D-EIRENE case (solid lines with points) and for the 5th iteration EDGE2D-EIRENE case (solid lines without points). Almost no difference between profiles of the two iterations can be seen in the figure, except for the 5th iteration T i profiles having higher values in the outer core and in the SOL. A significant difference between T i profiles belonging to different iterations, compared to a much smaller difference between T e profiles, is attributed to a much lower ion than electron dimensionless collisionality in the SOL due to much higher upstream T i (plots of the dimensionless collisionalities can be found in figure 3 of [4]). Figure 3 shows parallel profiles of T i , main ion density n i , and ion conductive power flux q cond i (which is =0 at targets in EDGE2D), for poloidal ring s01 just outside of the separatrix, for various EDGE2D-EIRENE iteration numbers. Heat flux limits push maximum upstream T i values higher by factor 1.20, from 164.3 eV for the 0th iteration to 197.4 eV for the 5th iteration, in response to falling χ i . In the zones of highest q cond i , multipliers used for χ i values in EDGE2D were ≈0.22 and ≈0.13 at the high field side (closer to the inner divertor) and at the low field side (closer to the outer divertor), respectively, in the 5th iteration case. Parallel ion conductive power fluxes to divertors are lower in the 5th iteration case, while n i at entrances to divertors are higher. Power balance in the SOL including target powers and volume radiation will be discussed in section 5.
For electron parameters, parallel profiles look qualitatively similar to the ions', but owing to lower dimensionless electron upstream collisionality, T e and q cond e variations among different iterations are smaller. In particular, maximum T e values rise only by factor 1.10, from 90.44 to 99.37 eV. In the zones of highest electron conductive power fluxes multipliers used for χ e values in EDGE2D were ≈0.60 and ≈0.45 at high and low field sides, respectively, in the 5th iteration case. These profiles are not shown here. The same profiles as shown in figure 3 for ions, but at the poloidal ring s11, are shown in figure 4. This ring corresponds to a position inside the SOL where the maximum, or close to maximum, target power loads are obtained, as will be shown below in this section. In difference to the ring s01, maximum T i gradients are formed further away from divertors in the main SOL. Maximum upstream T i values rise by factor 1.36, from 96.0 to 130.8 eV, from 0th to 5th iterations in response to falling χ i caused by heat flux limiting. In the zones of highest q cond i , multipliers used for χ i values in EDGE2D were ≈0.26 and ≈0.10 at the high and low field sides, respectively, in the 5th iteration case. Similar to the ring s01, q cond i to divertors are lower in the 5th iteration case, while n i in divertors are higher. Stronger upstream T i rise from the 0th to 5th iteration cases than at ring s01, as well as lower q cond i multipliers, are related to lower collisionalities for upstream parameters at this ring compared to the ring s01: despite lower upstream T i , the effect of a much lower plasma density at ring s11 compared to ring s01 results in lower dimensionless collisionality, as can be seen from figure 3 of [4] (note that ring s01 in this reference is denoted as position i = 6). At the same time, electron dimensionless collisionality, unlike ion's, is higher at ring s11 compared to ring s01. Correspondingly, electron kinetic effects at this ring are weaker than at ring s01. Such a disparity between ion and electron profiles can be explained by electron dimensionless collisionality weakly rising from the ring s01 to ring s11, while for ions it is sharply falling. This is related to flatter T i than T e profile across the SOL. Electron profiles, similar to those plotted in figure 4 for ions, are not shown here.
Qualitatively similar effects of kinetics on temperature profiles as those shown in figures 3 and 4, with upstream temperatures rising and downstream falling, were found in earlier kinetic calculations, for example in the particle in cell (PIC) code modelling, when comparing them to fluid code calculations, see figure 6 in [7].  Figure 5 shows T e , n e , and total, ion plus electron, conductive plus convective, parallel power flux profiles to divertor targets for EDGE2D-EIRENE iterations 0 and 5. The effect of the heat flux limiting upstream, partly compensated by the rise in upstream temperatures, has a net effect of the reduction of the total target power flux. The reduction is achieved mostly by the target T e reduction. The input power into the EDGE2D grid is the same for all iterations, hence, a drop in the target power must be compensated by the rise of volumetric power losses. The detailed analysis of the power balance in these EDGE2D-EIRENE cases will be discussed in section 5.

Challenges to the iterative KIPP-EDGE2D coupling
Success in testing the algorithm of the iterative KIPP-SOLPS coupling in [6] for kinetic electrons can be attributed to relatively smooth profiles of plasma parameters along the field line (a 1D SOLPS version was used in [6]), with not very large T e drops at targets and relatively small cell-to-cell variations of plasma parameters in this reference. These are conditions where KIPP, whose parallel free-streaming scheme is based on the high resolution 2nd order explicit scheme with monotonic central-difference flux limiter ('MC limiter'), can give reliable results. The scheme uses reconstruction-evolution-averaging (REA) finite volume algorithm described in [8]. In edge 2D fluid codes, at the same time, the situation is often quite different, with significant cell-to-cell parameter variations. This is also the case with EDGE2D-EIRENE solutions analysed here, with plasma parameters close to entrances to divertors exhibiting large cell-to-cell variations. The 'evolution' part of the REA algorithm describes advection of a piece wise interpolated distribution function f. The slope of segments inside of a velocity cell is equated to the spatial derivative of f calculated In extreme cases where f cell-to-cell variations vary by a large factor, probably an exponential interpolation onto cell faces instead of the linear one would have been more adequate, yielding much lower power fluxes than according to the standard REA algorithm, which is expected to over-estimate the fluxes. Figure 6 shows profiles of T e , n e , and total parallel power flux (ion, including impurities, plus electron, conductive plus convective) at ring s01 for the original, 0th iteration EDGE2D-EIRENE case. A steep T e drop, from 46.0 to 4.5 eV between neighbouring cells across the entrance to the outer divertor, is indicated on the figure. It continues for one cell further into the divertor (see figure 7 in the next section). Such a steep drop, also seen on the T i profile, as well as a steep rise in n e , violates conditions of applicability of the KIPP freestreaming procedure, making it strongly imprecise. One can also see unusual features on the total parallel power flux profile, which exhibits a significant drop, by factor 1.44, in the gradient of the total power flux approaching the outer divertor from the main plasma side: ∇q(s || ≈ 92 -95 m)| is lower than ∇q(s || ≈ 89 -92 m)| by this factor. Near the entrance to the inner divertor, the power flux profile is non-monotonic, with its absolute value increasing from the cell face in the main SOL to the adjacent cell face at the entrance to the divertor. The later feature is impossible to justify physically, and it demonstrates that the EDGE2D solver struggles to deliver correct profiles near divertor entrances.
Strongly non-rectangular cells, in particular, EDGE2D grid cells around the X-point (see figure 1), may have contributed to large T e , n e and total power flux variation seen in figure 6. They also present an additional numerical problem for KIPP whose equations employ the philosophy of a flux tube. When reading EDGE2D-EIRENE output data into KIPP, it is implicitly assumed that an infinitely narrow flux tube passes through centres of EDGE2D cells. In reality, EDGE2D solves equations in flux coordinates, so their reconstruction in the Cartesian coordinate system leads to inaccuracies in the mapping, resulting in wiggles in the profiles. The mapping inaccuracy depends on the local change of the grid size which is the largest near the X-point.
Large cell-to-cell plasma parameter variations are, however, not limited only to the ring s01, despite strongly shaped non-rectangular cells are only present on this ring. On ring s02, the steep T e drop starts from the 1st cell already in the outer divertor, where T e = 33.1 eV, and following cells (the first of which is indicated in figure 1), with T e for the next three cells being 22.40, 6.35 and 3.13 eV. For the next ring, s03, the steepest T e drop occurs between the 6th and 7th (indicated in figure 1) cells counting from the cell just before the entrance to the outer divertor on the main SOL side, from 9.46 to 4.70 eV. These features are caused by sharp ionization fronts taking place at different locations in the divertor at different rings.
Discussion about challenges facing KIPP cases under conditions of strong cell-to-cell plasma parameter variations will be continued in the next section.

Parallel power flux and temperature profiles along field lines
In this section parallel profiles of plasma temperature, density and parallel power fluxes for the two selected poloidal rings, s01 and s11, are presented. These are the same rings as were chosen for the analysis in the two earlier papers [3,4], in which, however, they were labelled as slice numbers i = 1 and 6, with slice numbers i = 1-6 representing only a fraction of the total 16 poloidal rings covering the EDGE2D grid in the SOL and divertor.
All profiles in this section are plotted along the spatial grid used by KIPP, which, as was pointed out above, for numerical reasons of smoothness of KIPP solutions, has evolved out of the EDGE2D grid by dividing each EDGE2D cell into two equal size cells and interpolating parameters onto the new, finer grid. The final EDGE2D iteration 5 case is used here for plotting conductive power fluxes extracted from EDGE2D, as well as for convective power fluxes, which depend on profiles of plasma temperatures, densities and parallel electron and main ion velocities obtained in EDGE2D cases. Conductive power fluxes from KIPP also refer to the last, 5th iteration KIPP run. Since the outer target receives higher power flux than the inner target, only zoomed up profiles near this target will be presented, as in [3,4]. in the main SOL, as a result of the iterative KIPP-EDGE2D coupling, can be seen. The only exception is around the stagnation point at maximum T e values where, due to different signs of EDGE2D and KIPP conductive power fluxes, KIPP/EDGE2D multipliers were not applied (they were specified as unities in EDGE2D) and the two conductive power fluxes strongly diverge, see section 2 for explanations. Note that, since the KIPP/EDGE2D multipliers were not applied to EDGE2D cells in the divertor, convergence between KIPP and EDGE2D cases was not reached in this region, and KIPP and EDGE2D power fluxes may differ. KIPP results in the divertor should be considered as representing kinetic solutions for given background (thermal) plasma parameter profiles.

Profiles along ring s01
Large spikes of q cond,KIPP e at divertor entrances should be considered as artefacts of the numerical scheme for the freestreaming used in KIPP, as was discussed in section 3. Note that T e and power flux oscillations in figure 7 are smoothed by the procedure of doubling the number of EDGE2D cells to obtain the KIPP grid, with a consequent interpolation of plasma parameters, see section 2. This explains, in particular, why q cond,EDGE2D e spikes near divertor entrances cover two cell faces instead of one. Figure 8 shows the same quantities as shown in figure 7, but at distances near the outer divertor target, to provide good spatial resolution. q cond,EDGE2D e values are extremely low in the divertor, since they scale with T 5/2 e ∇ || T e , with T e being very low near the target. At the same time, q cond,KIPP e is much higher, being only a factor 2 lower than q conv e . The relatively flat q cond,KIPP e profile when approaching the divertor is explained by this power flux being carried mostly by superthermal electrons originating from regions far upstream in the main SOL [3]. Since these electrons are much less collisional than thermal electrons, the heat flux they carry is strongly nonlocal and not being much affected by a steep T e drop right at the target. Lower q cond,KIPP e than q conv e at the target implies that kinetic effects cannot strongly increase electron power deposition at the target calculated by EDGE2D-EIRENE.   region around the stagnation point at maximum T i values, can be seen, for the same reason as explained for electron profiles above. A sharp increase in q conv,KIPP i and a simultaneous sharp drop in q cond,KIPP i at the entrance to the outer divertor can be explained by the sharp increase in the main ion velocity towards the divertor entrance, as follows from the EDGE2D-EIRENE solution. Figure 10 shows the same quantities as shown in figure 8, but at distances near the outer divertor target. Qualitatively, they show the same features as can be seen in figure 9 for electrons. In particular, the relatively flat q cond,KIPP i profile approaching the divertor is explained by this power flux being carried mostly by super-thermal ions originating from the main SOL [4]. As for electrons, target q cond,KIPP i values Figure 11. The same profiles as shown in figure 7, but for the poloidal ring s11. are somewhat lower than q conv i . Note that q cond,KIPP i values are comparable to q cond,KIPP e (see figure 8) in the divertor. This is despite the factor ∼ √ m e /m i of reduction of classical collisional conductive ion power fluxes compared to electron power fluxes caused by lower ion velocities (for equal ion and electron temperatures). Apparently, much lower ion dimensionless collisionality caused by much higher T i than T e upstream (T i,max = 160.3 eV, T e,max = 88.8 eV) leads to a very significant increase in the ion parallel power flux in the divertor due to much less collisional (compared to electrons) non-local transport of super-thermal ions.
Overall, results of this sub-section show that for the given JET experimental conditions kinetic effects of the parallel charged particle transport do not critically influence target power load along the ring s01 with large temperature variations and partially detached plasma near the targets. Figure 11 shows parallel profiles of the same quantities as shown in figure 7, but for ring s11. Compared to ring s01, the plasma is well attached to divertor targets, and the T e variation along field lines is quite moderate, from the maximum of 36.5 eV to 5.82 and 11.7 eV at inner and outer targets, respectively. As for the ring s01, good agreement between q cond,KIPP e and q cond,EDGE2D e can be seen in the main SOL, except for the region surrounding the stagnation point with maximum T e . The convective power flux q conv e plays a secondary role in the total electron power flux, which is almost always the case for electrons, except for very cold (low T e ) and dense (high n e ) regions in the divertor. The negative value of q conv e almost along the entire ring reflects the plasma (ion) flow from the outer to inner target in the EDGE2D-EIRENE solution. Figure 12 shows zoomed up parallel profiles of the same quantities as shown in figure 11, for the region near the outer target. A relatively modest T e variation along field lines,   Figure 13 shows parallel profiles of the same quantities as shown in figure 9, but for ring s11. As with all other conductive power profiles in the main SOL presented above, good agreement between q cond,KIPP i and q cond,EDGE2D i in the main SOL can be seen, except for the region around the stagnation point with the highest T i . Ion temperature along field lines is substantially higher than electron, varying between T i,max = 130.8 eV and 19.0 and 8.54 eV near inner and outer targets, respectively. A strong decoupling between T i and T e is caused by low dimensionless collisionalities of both species, including in the divertor region. An unusual situation with the substantially lower T i at the outer than at the inner target, following from the EDGE2D-EIRENE solution, must be caused by the reversal of q conv i inside of the outer divertor, transferring thermal energy 3/2n i T i to kinetic energy n i m i V 2 ||i /2. Equipartition energy exchange between ions and electrons, on the other hand, cannot be responsible for lower T i at the outer target due to T e > T i near the outer target, as can be seen from figures 12 and 14 (see below). Note that in the main SOL the convective power flux q conv i plays a much larger role in the total ion power flux than for electrons, as expected from relatively low ion parallel velocities compared to electron velocities (factor ∼ √ m e /m i , for comparable T e and T i ). Figure 14 shows zoomed up parallel profiles of the same quantities as shown in figure 13, for the region near the outer target. The dominant role of ion convection in the total parallel target ion power flux can be seen. The much higher q cond,KIPP i compared to q cond,EDGE2D i must be caused by a much lower ion dimensionless collisionality at this ring, leading to a large contribution of ion non-local kinetic effects to the conductive power flux. This is related both to a very high T i,max = 130.8 eV (compared to T e,max = 36.5 eV) and to a very low upstream dimensionless collisionality. According to the ion and electron upstream (at the OMP position) dimensional collisionality profiles, which can be found in figure 3 of [4], the ion dimensionless collisionality on ring s11 is lower than the electron dimensionless collisionality by factor ≈10.

Profiles along ring s11
Despite very strong heat flux enhancement (a very large contribution of ion non-local kinetic transport effects leading to a much higher parallel ion conductive power flux compared to that calculated in EDGE2D according to a strongly collisional fluid model) the total ion power flux to the outer target, with the inclusion of the convective power flux contribution, appears to be not too much affected by ion kinetic effects, as was the case with the ion target power flux for ring s01 and for electron power fluxes for both rings s01 and s11.

Power balance in KIPP-EDGE2D solutions
Contributions of electron and main ion conductive and convective parallel power fluxes at the outer target, with the exception of the ion conductive power flux from EDGE2D, q cond,EDGE2D i , which is zero, for the 5th iteration of both EDGE2D-EIRENE and KIPP runs, are plotted in figure 15. This figure can be compared with figure 6 of [4], which shows some of the fluxes presented here for selected rings for what is here referred to as the 1st KIPP iteration (following the 0th EDGE2D-EIRENE iteration). It has to be pointed out, however, that in the legend of figure 6  In the region of highest target power loads the largest power flux is kinetic conductive electron power flux q cond,KIPP e . It exceeds the fluid power flux q cond,EDGE2D e by factor ≈2 owing to the contribution from the non-local transport of super-thermal electrons. Indicating 'heat flux enhancement' due to the contribution from the non-local transport of super-thermal electrons. Note however, that q cond,EDGE2D e is lower than the sum of other power flux contributions.
In the outer SOL/divertor region the largest power flux comes from ion convection, q conv i . The dominance of this term, and in particular, its access over q conv e , is caused by much higher T i than T e in the whole SOL and divertor region owing to much flatter upstream T i profiles. This region is also characterized by small T e,i variations along field lines, leading in particular to fairly close values of electron conductive power fluxes calculated by KIPP, q cond,KIPP e , and EDGE2D, q cond,EDGE2D e .
In none of the above regions ion conductive power flux, even with inclusion of effects of non-local transport of superthermal ions, q cond,KIPP i , plays any significant role. As was discussed in section 2, kinetic effects in the main SOL (leading to 'heat flux limiting') reduce the total target power flux, as was demonstrated by comparing total power fluxes extracted from 0th and 5th iterations of EDGE2D-EIRENE cases shown in figure 5, with the exception of power fluxes at the outer target on rings s12-s16 (the last five points in figure 5) connected to the outer SOL. At these outer rings small T e,i variations along field lines and, consequently, small contributions of kinetic effects to power fluxes follow from KIPP. Lower total target power load, by factor ≈1.4, in the 5th iteration EDGE2D-EIRENE case, 0.771 MW, compared to the 0th iteration case, 1.09 MW, is confirmed by the EDGE2D output. Lower power load in the 5th iteration EDGE2D-EIRENE case is compensated by higher volumetric power losses, mainly from higher electron power loss on radiation and higher ion power loss on charge exchange. These increased volumetric power losses must be caused by higher plasma density in zones of maximum radiation. Higher n i in the 5th iteration EDGE2D-EIRENE case can be seen around the X-point for ring s01 (figure 3), and closer to the target for ring s11 (figure 4), with higher plasma densities in these regions in turn being caused by lower ion and electron temperatures.
In the iterative KIPP-EDGE2D-EIRENE modelling described here non-Maxwellian features of the electron distribution function in the calculation of rate constants for atomic physics processes contributing to particle and power sources were not taken into account, despite a significant impact of high energy tail electrons on the rates in the divertor was reported (see e.g. [7,[9][10][11]). The main trends in adding kinetic atomic rates to fluid code plasma solutions in high recycling plasmas are summarized in [12], where 1% population of high energy super-thermal electrons was artificially added to SOLPS (a 2D edge fluid code) modelling into detached divertor plasma. It was shown that the presence of high energy electrons in the low T e ≈ 1 eV thermal population, high density plasma caused an enhanced ionization rate; then, since the recombination rate was essentially unchanged, ions recombined. Thus, a new 'volume recycling' zone was established, with enhanced energy losses, and ultimately, less ions reached the target, assisting the detachment process.
In KIPP calculations presented here energies of superthermal electrons on field lines with detachment are extremely high, ∼a few 100's of eV, since they mostly originate from the regions of the main SOL with highest T e , with their energies being ∼6T e , as was estimated in [3]. At the same time, the rate coefficients for atomic and molecular hydrogen ionization have a maximum at or below T e ≈ 100 eV (see e.g. [13], figure 1.25). The cross-section for nitrogen ionization also peaks at T e ≈ 100 eV ( [14], figure 3). Hence, no particularly strong impact of ionization/excitation processes on super-thermal bump-on-tail electrons is expected. And in any case, such processes are expected to increase the ionization rate and detachment in the divertor. The results obtained in the present study, pointing to only a moderate increase in target power fluxes caused by kinetic effects, can therefore be regarded as conservative in the sense that the proper description of the atomic rates may only increase the volumetric power loss and the degree of detachment, thereby reducing the target power load.

Summary
Owing to large T e and T i variations along field lines close to the separatrix and partial detachment at divertor targets, the JET H-mode radiative divertor conditions analysed in this paper may to some extent be considered reactor relevant. Further away from the separatrix conditions change, from strongly attached, with high power flux conducted to the target, to low recycling conditions in the outer SOL, with much lower upstream T e,i and relatively flat temperature profiles along field lines.
The main progress achieved in this work compared to earlier studies [3,4] is in the iterative coupling between KIPP and EDGE2D, and with KIPP cases running for both ions and electrons simultaneously, resulting in the converged KIPP-EDGE2D-EIRENE solution, Despite the failure of the iterative coupling scheme between KIPP and EDGE2D in the divertor region owing to very high T e , T i and n e cell-to-cell variations in EDGE2D-EIRENE solutions, some conclusions, based on the successful performance of the iterative scheme in the main SOL, can be made. The results of this work, combined with those of earlier studies in [3,4], for JET conditions analysed here can be summarized as follows: • On field lines with large T e and T i variations in the main SOL, kinetic effects of non-local parallel transport of superthermal electrons and ions reduce conductive power flows, leading to 'heat flux limiters' with reduced parallel conductive power flux coefficients compared to fluid code solutions. This is a well-known result of kinetic calculations (see e.g. [3] and references therein). As a result, maximum upstream temperatures, T e,i,max , increase, but this does not restore original (following from fluid calculations) parallel power fluxes. This in turn leads to the reduction of power fluxes to the divertor, the reduction of T e,i and the increase of plasma density and volumetric power losses due to plasma interactions with neutrals and impurities in the divertor. • In the near SOL, on field lines just outside of the separatrix, with partial detachment near strike points, owing to very low target T e,i dominant contributions to the outer target power flux come from electron and ion convection. High heat flux enhancement factors for conductive power fluxes caused by kinetic effects of parallel propagation of super-thermal electrons and ions do not make these fluxes large enough to compete with convective power fluxes, as classical collisional conductive power fluxes are almost negligible due to very low T e,i in the divertor. The electron conductive power flux calculated by KIPP may, however, be close to the electron convective power flux, which is attributed to the power flux carried by a tiny minority of super-thermal electrons originating from locations in the main SOL with highest T e . • In the middle of the SOL, on field lines with attached divertor conditions and with highest target power loads, the dominant contribution to the outer target power load comes from electron conduction, which, due to kinetic effects, is higher by factor ∼2 compared to that following from EDGE2D-EIRENE, which, in turn, is close to convective power fluxes of both electrons and ions. The contribution from the ion conductive power flux, even despite its strong enhancement by kinetic effects, is insignificant. • In the outer SOL the dominant contribution to the outer target power load is from the ion convective power flux, which is factor ∼2 higher than the electron convective power flux owing to much higher T i than T e . Conductive and convective electron power fluxes are close to each other. Ion conductive power flux is negligible.
A very important conclusion from results of this work, which is in line with earlier studies [3,4], is that under inter-ELM H-mode JET radiative divertor conditions analysed, kinetic effects of parallel charged particle propagation do not significantly increase the total outer target power load near the strike point. The kinetic effects are unlikely to cause a power flux 'burnthrough' leading to a substantial increase in the target power load at locations where the fluid code (EDGE2D-EIRENE) predicts partial divertor detachment. On the contrary, the impact of a minority of super-thermal electrons on atomic rate coefficients can only increase the volumetric power loss in the divertor, leading to stronger detachment, as pointed out in section 4.
For field lines with attached plasma and a relatively high power flux to the divertor target (in the middle of the SOL), kinetic effects of parallel electron transport may raise the conductive electron power flux by factor ∼2, which then becomes the dominant contribution to the total target power flux. But even there, the total target power flux is increased by only factor ∼1.5 compared to the total target power flux predicted by EDGE2D-EIRENE calculations, when contributions of convective electron and ion power fluxes are taken into account.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.