EDGE2D-EIRENE simulations of the influence of isotope effects and anomalous transport coefficients on near scrape-off layer radial electric field

EDGE2D-EIRENE (the ‘code’) simulations show that radial electric field, Er, in the near scrape-off layer (SOL) of tokamaks can have large variations leading to a strong local E × B shear greatly exceeding that in the core region. This was pointed out in simulations of JET plasmas with varying divertor geometry, where the magnetic configuration with larger predicted near SOL Er was found to have lower H-mode power threshold, suggesting that turbulence suppression in the SOL by local E × B shear can be a player in the L–H transition physics (Delabie et al 2015 42nd EPS Conf. on Plasma Physics (Lisbon, Portugal, 22–26 June 2015) paper O3.113 (http://ocs.ciemat.es/EPS2015PAP/pdf/O3.113.pdf), Chankin et al 2017 Nucl. Mater. Energy 12 273). Further code modeling of JET plasmas by changing hydrogen isotopes (H–D–T) showed that the magnitude of the near SOL Er is lower in H cases in which the H-mode threshold power is higher (Chankin et al 2017 Plasma Phys. Control. Fusion 59 045012). From the experiment it is also known that hydrogen plasmas have poorer particle and energy confinement than deuterium plasmas, consistent with the code simulation results showing larger particle diffusion coefficients at the plasma edge, including SOL, in hydrogen plasmas (Maggi et al 2018 Plasma Phys. Control. Fusion 60 014045). All these experimental observations and code results support the hypothesis that the near SOL E × B shear can have an impact on the plasma confinement. The present work analyzes neutral ionization patterns of JET plasmas with different hydrogen isotopes in L-mode cases with fixed input power and gas puffing rate, and its impact on target electron temperature, Te, and SOL Er. The possibility of a self-feeding mechanism for the increase in the SOL Er via the interplay between poloidal E × B drift and target Te is discussed. It is also shown that reducing anomalous turbulent transport coefficients, particle diffusion and electron and ion heat conductivities, leads to higher peak target Te and larger Er, suggesting the possibility of a positive feedback loop, under an implicitly made assumption that the E × B shear in the SOL is capable of suppressing turbulence.


Introduction
Scrape-off layer (SOL) and divertor make a direct impact on the core region of the plasma via neutral penetration and ionization. While ionization of neutrals is necessary for achieving the desired plasma density in the core, it was known long ago that high large neutral density can lead to the degradation of plasma energy confinement. As was stated in the work by Wagner et al in 1985, in reference to experimental results of ASDEX and PDX tokamaks in achieving transition from low to high energy plasma confinement regimes (L-H transition) [1]: 'The detrimental effect of enhanced recycling and gas puffing has been experimentally observed. Both PDX and ASDEX have observed that gas fuelling from the divertor dome instead of the main plasma chamber facilitates the H transition. Excessive gas puffing can quench the H-mode of PDX yielding low confinement'. Other (than gas puffing) methods of plasma fuelling, avoiding excessive neutral density at the plasma edge, such as neutral beam injection and deep pellet fuelling, are found to result in desired core density levels without causing confinement degradation. As was recently pointed out by Valovic et al [2]: 'A reasonable working hypothesis is that the optimum way to fuel the plasma across the L-H threshold is to increase the plasma density inside the separatrix while minimizing the number of particles deposited in the scrape-off layer'.
Recent advance in measurement techniques of plasma edge parameters showed that even when neutrals' ionization does not significantly alter plasma parameter profiles in the core, changing divertor geometry can still make a strong impact on the H-mode power threshold, P LH , (factor two variation in the case of JET plasmas in the ITER-like Wall (ILW) environment, with W divertor and Be main chamber [3][4][5]). The divertor geometry effect on confinement was also seen as the effect of the X-point height (vertical distance between the X-point and the material surface) in JET [6,7], MAST [8], DIII-D [9], Alcator C-mod, where it was referred to as the effect of the outer divertor leg length [10], and COMPASS [11]. All these experiments suggest that there may be other, more subtle mechanisms for the influence of the SOL and divertor on the core plasma, rather than the direct effect of the edge density rise due to neutral ionization. Here it is worth noting that, while the main trends in the parameter dependence of the H-mode power threshold scaling are usually related to core plasma (input power, average density) and engineering parameters (machine size, magnetic field), it is well known that the L-H transition physics is an edge physics process. In the multi-machine Martin's scaling for the H-mode power threshold [12] there is a large, factor four scatter in experimental data points, apparently due to unaccounted factors/parameters. As pointed out above, only the variation of divertor geometry in JET along causes the factor two variation in P LH .
In search of unaccounted parameters which may explain a large scatter of experimental points for the H-mode power threshold, EDGE2D-EIRENE code [13][14][15] was used to simulate two JET discharges with different divertor configurations but very similar plasma parameter profiles in the core. It has to be noted that, while EDGE2D-EIRENE is often referred to as a '2D fluid edge code', in reality it is a code package which contains coupled fluid EDGE2D code for plasma behavior coupled to the Monte Carlo code EIRENE for neutrals.
In one of the simulated discharges the outer strike point was on the horizontal divertor target (HT configuration), while in the other-on the vertical divertor target (VT configuration). Inner strike points were both on vertical targets. Experimentally, in the HT configuration the H-mode power threshold was found to occur at about twice less input power than in the VT configuration. In EDGE2D-EIRENE simulations, the largest difference between the two cases was seen in the E r profiles at the outer midplane (OMP), originating due to very different target T e profiles, see figure 4 in [16]. In this figure one can see a very large E r spike in the near SOL in the HT configuration, and it was hypothesized that the strongly localized E×B shear in the SOL could suppress turbulence resulting in lower H-mode threshold power. Such E r spikes were not observed in the experiments. Measurements of this quantity are currently planned in JET and ASDEX Upgrade.
A substantial difference between near SOL E r spikes was found in EDGE2D-EIRENE cases with different hydrogen isotopes, with largest spikes for heavier isotopes (the highest -in tritium). This correlates with experimental observations of the lower H-mode power threshold in plasmas with larger average isotope mass [17]. Isotope experiments and modeling will be discussed in more detail in the next section.
Finally, a very different approach to a possible influence of neutral ionization, including the impact of the neutrals' isotope mass and the influence of the target material on their energies, compared to the one considered in this paper, should be mentioned. Part of the radial electric field known to be attributed to the non-zero angular momentum can be reduced by charge-exchange (CX) neutrals. In the study [18] it is demonstrated that the neutral isotope mass influences the momentum loss, and, via it, E r , in a way which is qualitatively consistent with experimental observations of better energy confinement in plasmas with heavier hydrogen isotopes. In addition, higher Z wall materials result in higher neutral energy reflection coefficients from the wall, deepening penetration of neutrals into the plasma. The notion 'fuelling fuels turbulence' is introduced to describe detrimental effects of neutral CX processors on the confinement [18].
In this paper section 2 contains a brief overview of isotope effects on plasma confinement, relevant to the present study. Setup of EDGE2D-EIRENE modeling cases is discussed in section 3. Results of EDGE2D-EIRENE cases are discussed in section 4, including upstream E r and target parameter profiles (section 4.1), neutral ionization pattern (section 4.2), drift effects (section 4.3), and the influence of varying anomalous transport coefficients on near SOL E r (section 4.4). Conclusions from the work are drawn in section 5.

Isotope effects on plasma confinement
As was pointed out earlier in this paper, the H-mode power threshold is found to be lower in plasmas with heavier hydrogen isotopes. The plasma confinement is typically also worse in H compared to D plasmas (most of the experiments were done in hydrogen and deuterium plasmas). Since modeling results presented in this paper center on effects of different hydrogen isotopes it is worth giving a brief overview of experimental findings, limited mostly to ASDEX Upgrade and JET L-mode plasmas, and simulations (covering also H-modes, as these results are scarce) with 2D edge codes SOLPS and EDGE2D. Both codes use EIRENE for the description of neutrals behavior.
One of the first attempts to simulate ASDEX Upgrade experiments, focusing on the difference between plasma edge parameters in large number of H and D plasmas, was made with a variant of SOLPS, B2.5-I, adapted for automatic changes in transport coefficients to match experimental density and temperature profiles at the plasma edge, including SOL [19]. It was found in this study that in H-mode discharges heat conductivities were lower in D plasmas. A thorough analysis of two ASDEX Upgrade H-mode discharges in H and D was carried out in [20] using SOLPS simulations. Perpendicular anomalous diffusion coefficient, D ⊥ , and electron heat conductivity, χ e , were established to be lower in the D plasma in both outer core and SOL regions, while ion heat conductivities, χ i , were the same in H and D plasmas. Minima in transport coefficients were reached at a distance~1 cm inside of the separatrix (at the OMP position). Recently, EDGE2D-EIRENE interpretative simulations were used to compare JET-ILW NBI heated L-modes in H and D in the near separatrix region [21]. Again, larger D ⊥ was established in the H discharge, less information was available on heat conductivities χ ei due to poorer quality of experimental temperature profiles.
In recent ASDEX Upgrade L-mode experiments with ECRH heating and T e >T i it was established that the isotope dependence of energy confinement is a consequence of the mass dependence in the collisional equipartition energy exchange between electrons and ions, as confinement in hydrogen-compared to deuterium-plasmas is accompanied by a larger ion heat flux due to the enhanced electron to ion heat transfer, while the electron transport channel is unaffected by the difference in ion mass [22]. This explanation however is not applicable to JET-ILW NBI heated L-modes, where ion and electron temperatures are close to each other and the equipartition energy exchange cannot be computed.
Recent EDGE2D-EIRENE simulations of JET ELMy H-mode discharges in H and D with similar stored energies and line-average densities also showed stronger transport barriers for particle and energy fluxes, with transport coefficients being lower in the D discharge. In both cases the minima in transport coefficients were achieved just inside of the separatrix. These discharges however were performed in somewhat different magnetic equilibria.
Summarizing, we conclude that available 2D edge code simulations support experimental finding for poorer energy and particle confinement in H compared to D plasmas which is most likely caused by the differences in confinement of the edge plasma. They also show that anomalous transport coefficients are larger in H plasmas not only in the core, but also in the near separatrix and SOL regions. Regarding a possible impact of the near SOL E rr on the plasma confinement described in section 1, there is however no (mostly indirect) evidence that such a link really exists, unless the plasma is close to the L-H transition boundary, where the E×B shear in the near SOL could potentially accelerate the transition. At the same time, in plasmas with high input power, in particular in H-modes, matching stored energies in H and D requires substantially higher input power in H, leading to larger power flux through the separatrix, higher strike point T e in the divertor and larger near SOL E r , indicating the absence of the correlation between the near SOL E r and confinement. The isotope dependence of the plasma confinement should therefore be attributed to the pedestal physics.

Setup of EDGE2D-EIRENE cases
Setup of EDGE2D-EIRENE cases was similar to that described in [23], including the grid which was built based on the magnetic equilibrium of JET pulse #81883 with the outer strike point on the horizontal target (the grid is shown in figures 1 and 2 of [23]). All cases were run with drifts and parallel currents. A self-consistent neoclassical model for E r was implemented in the core which impeded surface averaged radial currents [24]. Only pure hydrogen isotope cases (in H, D and T) were run, without impurities. The material surfaces were assumed as in the ITER-like Wall (ILW). The EIRENE version with Kotov-2008 model was used to describe neutral behavior. Spatially constant anomalous transport coefficients were specified: diffusion coefficient D ⊥ =1 m 2 s −1 and electron and ion heat conductivities χ e,I =2 m 2 s −1 across the whole grid, unless otherwise stated. The input power into the grid was set at 3.0 MW, split equally between ion and electron channels.
In the earlier modeling described in [23] the plasma density was controlled by a combination of gas puff and wall recycling ('puff+recycling' option in EDGE2D-EIRENE), aiming at maintaining a specified electron density at the OMP position of the separatrix, n e,sep . The 'puff+recycling' option uses puffing if the desired n e,sep is above that specified by the modeler, whereas the reduction of the target recycling is used to reduce n e,sep if no puffing is applied. In the steady state the amount of gas puffed into the plasma has to be matched by the amount of gas pumped out by artificial pump surfaces placed in the divertor. In contrast with the earlier modeling, in the modeling described here the fixed gas puff from the top of the machine, at 4.2×10 21 s −1 , was used. This level of the gas puff corresponded to the case in deuterium described in [23]. The change from fixing n e,sep to constant gas puff for all isotopes was done in order to get an insight into physics governing particle confinement in plasmas with different hydrogen isotopes. Since in the experiment H discharges are found to require larger gas puff to reach the same line-average density as in D discharges, indicating poorer particle confinement in H, the question arises whether for the same gas puff in H as in D the near SOL E r is smaller in H. In such a case the hypothesis that the E×B shear in the SOL can suppress turbulence leading to better particle confinement would be strengthened. In earlier isotope scans, described in [23], gas puff levels had to be increased going from T to H in order to maintain the same n e,sep . This caused somewhat higher radiation power losses for lighter isotope plasmas, contributing to lower peak target T e and near SOL E r . The present setup is therefore more challenging for obtaining higher near SOL E r in heavier isotope plasmas.
Finally, the newest version of EDGE2D-EIRENE was used in this study, which however caused only very minor changes in output profiles of plasma parameters compared to the earlier version used in cases described in [23].

EDGE2D-EIRENE cases and analysis
There were two main motivations for running EDGE2D-EIRENE cases in this study. The first was to establish whether in heavier hydrogen isotopes, for the same level of the gas puff, the spike of the near SOL E r is larger, and if so, to understand the physics of this phenomenon. The second motivation was to test the effect of varying transport coefficients on the near SOL E r . If the near SOL E r is indeed larger in heavier isotopes, one would then expect stronger turbulence suppression in the SOL by the E×B shear and an improvement in particle and energy confinement. This would then justify using lower transport coefficients in cases with heavier isotope plasmas. At the same time, turbulence suppression leading to lower transport coefficients would result in shorter plasma density and temperature decay lengths thereby increasing E r . Hence a positive feedback loop between the near SOL E r and target T e , profile which causes it, and turbulence suppression by the E×B shear can be created. Such a positive feedback loop may lead to more substantial differences between plasmas in H, D and T compared to those caused by direct consequences of the variation in neutral penetration lengths into the plasma. In this study the suppression of turbulence in the SOL by the increase in the local E×B shear is considered as a hypothesis, since no modeling supporting this is presented.

OMP E r and target profiles in isotope scan cases
In this section the output from EDGE2D-EIRENE isotope scan cases with exactly the same setup are described. In particular, H, D and T cases were run with the same transport coefficients and the same gas puff. Figure 1 shows OMP profiles of n e , T e and T i versus distance from the separatrix. The largest difference between H, D and T cases is in density profiles, especially near the separatrix position and in the SOL. Separatrix density n e,sep is higher in the T case compared to the H case by factor 1.187. It is not clear however why n e,sep in the D case is much closer to n e,sep in the H case than to the T case, the reason for this unexpected result will be investigated. At the same time, charged particle fluxes through the separatrix, mostly attributed to the anomalous diffusive particle flux density −D ⊥ ∇n, equal to 1.89 × 10 22 , 1.29 × 10 22 and 1.06 × 10 22 s −1 for H, D and T cases, respectively, exhibit the expected tendency, closely following the inverse square root dependence on the isotope mass. This can be explained by deeper penetration of lighter neutrals into the plasma and their stronger ionization source in the core. Since most of the ionization takes place within the simulation domain (EDGE2D grid), in the steady state conditions, achieved in each of the cases, ionization source in the core must be compensated by the plasma particle flux across the separatrix from the core into the SOL, which is largest in the H case.
OMP E r profiles for all three cases are shown in figure 2. There is a very little difference between these profiles in the core, and E r is almost constant for all cases there. The similarity of E r profiles is expected since the neoclassical E r doesn't depend on the ion mass, see the discussion below. The only significant feature of the E r profiles, capable of creating large E×B shear, is in the SOL, and in particular, positive E r spikes in the near SOL. The spike is largest in the T case and lowest in the H case. However, positives spikes in the near SOL for D and T are very close to each other, while the negative spike further out in the SOL is largest for H.
Profiles of E r in the core region are rather flat. This is due to flatness of ion pressure and T i profiles which determine E r in the neoclassical theory, provided there is no toroidal momentum (which was assumed to be zero at the inner core boundary for all cases presented in this paper). Plasma transport equations in EDGE2D assume strong collisionality Pfirsch-Schlüter regime for both SOL and core plasma conditions. In the core, however, this assumption is at least marginal, and is often violated, with the collisionality regime becoming Plateau. The coefficient K which accounts for the contribution of the T i gradient to the neoclassical E r , given by K∇T i /e, depends on the collisionality regime, being the largest in the Pfirsch-Schlüter regime, 2.1, smaller in the Plateau regime, 1.5, and even smaller in the Banana regime, −0.17 (see e.g. [25]). This is in addition to the ion pressure gradient contribution ∇ pi /n e . If the EDGE2D model was more realistic, by covering also the Plateau regime, it would yield steeper E r profiles, with lower (larger negative) E r just inside of the separatrix. Another reason for flatness of EDGE2D E r profiles is the use of spatially constant anomalous transport coefficients in the runs described here. In the experiment these coefficients are usually lower just inside of the separatrix leading to higher E r . Such changes in the code would have made its E r profiles qualitatively more consistent with experimental observations of the E r 'well' just inside of the separatrix (see e.g. ASDEX Upgrade results in [26]).
Target profiles of n e , T e and T i for all three cases are shown in figures 3(a)-(c), respectively. Target densities are lowest in the H case, which can be explained by higher velocity of H neutrals and lower ionization source in the SOL (see analysis of ionization sources in the next sub-section). This also makes its impact on n e,sep , which is lowest in the H case, as discussed above. Peak T e at the outer target (OT) is highest in the D case. This does not contradict to the fact that the highest peak E r is seen in the T case, since it is the radial T e gradient, rather than the T e itself, which gives the largest contribution to the upstream E r , as follows from the explanation of the origin of E r based on the effect of the Debye sheath potential drop~3T e /e [27], under conditions of constant target potential (the target in these EDGE2D-EIRENE cases was assumed to be made of tungsten, hence, the target potential was radially constant). Since electron temperatures at the inner target (IT) are very low near the strike point, they make a negligible contribution to the formation of the upstream E r . At the same time, farther into the SOL, at distances>1 cm from the inner strike point magnitudes of T e spikes for different isotopes are in reverse order compared to positions closer to the (outer) strike point. These, inner target T e spikes have opposite radial dependence (T e is rising away from the strike point) and result in negative E r spikes at the OMP, which are largest for lighter isotopes. Target T i profiles are qualitatively similar to the target T e profiles. Judging from the values of E r spikes, positive spikes seem to be more important than negative ones and are capable of generating higher E×B shear.

Neutral ionization pattern in isotope scan cases
The total ionization source for each poloidal ring, summed up over all cells, from outer to inner target for SOL rings, and along the whole poloidal transit for core cells, versus OMP cell center positions relative to the separatrix position, is plotted in figure 4(a), and its expansion around the separatrix position is plotted in figure 4(b). Calculation of these quantities requires multiplication of ionization density (in m −3 s −1 ) by the cell volume, cell by cell, and summing up the results over all cells. Neutral ionization occurs mostly in the divertor. Ionization source in the core is largest in the H case, as expected, and consistent with n e profiles shown in figure 1.
The characteristic penetration length of neutrals into the core is impossible to calculate using quantities described just above owing to the 2D nature of the ionization pattern, with maxima of the total ionization reached near the X-point position due to the large flux expansion there. Deeper in the core, a long tail of neutral ionization can be seen in figure 4(a), which must be caused by more energetic CX neutrals. The initial source of neutrals in the system, in addition to the gas puff, is mostly recombination at the target, with volume recombination being much smaller. For a D case, the volume recombination rate is equal to 1.63×10 21 s −1 , which has to be compared with the total ionization rate of 2.81×10 23 s −1 . At the same time, local ionization rates along the OMP position obey the expected isotope mass dependence: e-folding decay lengths across the core region are 4.6, 3.1 and 2.4 cm for H, D and T cases respectively. Figures 5(a), (b) show poloidal distributions of ionization sources against cell row numbers for the first ring in the SOL (ring S01, according to the EDGE2D nomenclature, figure 5(a)) and for the last ring in the core, closest to the separatrix (ring C01, according to the EDGE2D nomenclature, figure 5(b)). Row numbers give numbers of cells in the poloidal direction. For SOL rings numbering starts from the OT, while for core rings it starts from the bottom of the grid. In either case the poloidal cell numbering adopted in EDGE2D is counter-clockwise. Profiles shown in figures 5(a), (b) give total ionization sources for each cell, including the cell volume. Maxima of ionization sources seen in figure 5(a), especially large in the H case, correspond to cell positions close to the X-point. The same is true for profiles shown in figure 5(b). Note that the number of cells is larger for ring S01, since this ring includes cells in the divertor in difference to the ring C01 which covers the core. Poloidal profiles of the   same quantity for rings S02 and C02, further away from the separatrix in the SOL and core, respectively, show qualitatively similar but less pronounced patterns.
The reason for large contributions from cells near the X-point can be seen from figure 6 which shows the expanded view of the EDGE2D grid in the divertor and the bottom of the core region. Due to large poloidal flux expansion around the X-point, EDGE2D grid cells around this position are larger than other cells. Consequently, the ionization source is largest in these cells. Among the three isotopes, hydrogen atoms, having the highest velocity, have largest probability to reach these cells, while for D and T larger proportion of neutrals gets ionized before reaching these cells. This explains why the H case has stronger ionization in the near SOL rings, adjacent to the separatrix, and weaker ionization in the far SOL, which gives rise to flatter T e profiles along OT and lower E r spikes at OMP.

Drift effects in isotope scan cases
Understanding certain features of target profiles shown in previous sub-sections require consideration of drifts effects, especially the effect of the poloidal E×B drift caused by radial electric field and directed from the IT to OT for positive E r , which is correct for the 'forward' toroidal field direction with the ion ∇B drift down, towards the X-point, which was the case in all EDGE2D-EIRENE cases discussed in this paper. In the presence of drifts, some commonly used assumptions, such as conservation of total, ion plus electron pressure along field lines in the absence of parallel momentum losses, or the law for parallel heat flux, where T e is upstream (e.g. OMP) electron temperature and L || -parallel distance between the upstream position and the target [27], do not apply. This is because poloidal E×B drift carries both particle and (convective) power fluxes towards the OT (for positive E r and 'forward' toroidal field direction). Target profiles in the hydrogen scan cases were strongly affected by poloidal E×B drift. Figure 7 shows an expanded view, zoomed around the outer strike point position, of the E r profile at the OMP and n e and T e profiles at OT of the isotope scan cases versus distance from the separatrix position mapped to the OMP. There is one missing point on the OMP E r profile, at the first SOL ring S01. E r points at this location are also absent in figure 2, but there it is less noticeable. The reason why these points aren't plotted is related to the procedure of the calculation of E r values in the SOL and in the core. In the SOL, electric potential is the primary quantity. It is possible to calculate it since there is a reference point for the potential, which is assumed to be zero at (conducting, as in these cases) divertor targets. E r is then calculated by subtracting electric potentials from the neighboring rings. In contrast, in the core the primary quantity is E r , which has values that ensure ambipolarity of drift-related surface averaged radial (normal to flux surfaces) plasma fluxes. Since core and SOL regions are topologically separated from each other, it is impossible to connect electric potentials in them, unless one uses some model for electric currents through the separatrix. Such a model is not used in EDGE2D, instead, it is assumed that electric potentials in the last (outermost) core ring C01 and the first (innermost) SOL ring S01 are equal to each other, hence E r across the separatrix is assumed zero. Since the calculation of E r for the SOL ring S01 requires subtraction  between electric potentials in rings C01 and S02, taking into account a rather arbitrary connection between potentials in the core and SOL, such a procedure does not give physically sensible E r values. For this reason E r values for the ring S01 are not plotted. It has to be stressed that E r profiles in the SOL shown in figure 2 for positive distances from the separatrix position are due entirely to SOL physics.
The primary quantity which determines relative changes in E r and electron density n e for the three isotopes shown in figure 7 is electron temperature T e . It shows a much flatter profile for H than for D and T for the first few rings just outside of the separatrix. Large E r differences between these cases are mostly the consequence of the differences in target T e profiles. Since positive E r values in the SOL result in poloidal plasma fluxes towards the OT (in these cases, for the give toroidal field direction), they cause the plasma compression at OT. This leads to larger n e at OT. The E×B compression also causes the T e rise at OT, which is explained just below. Owing to the coefficient 5/2 (instead of 3/2) in the expression for the convective power flux, G / T 5 2 , where G is particle flux and T is temperature, which equally applies to parallel and cross-field fluxes, plasma compression at the OT by the poloidal E×B drift causes increase of both plasma density and temperature. It can be easily shown that the effect of such a compression gives a fairly high rate of the temperature increase, equal to 2/3's of that for density: This, in addition to the tendency for lower target n e in the H case discussed earlier, further increases target T e in the first two poloidal rings in D and T cases compared to the H case. There may therefore exist a positive feedback loop between the near SOL E r and target T e , where poloidal E×B drift increases target T e , which, in turn, increases E r in the SOL. It has to be noted that T e spikes in the near SOL typically result in a spatially variable E r , especially if one considers a strong E r jump across the separatrix. Hence, these features should lead to the formation of the shear of the poloidal E×B velocity. At the same time, the increase of the target density as a result of the plasma compression at the target due to the E×B drift leads to larger ionization and radiation near the target, which may limit the T e rise and thereby limit the feedback effect. The impact of the poloidal E×B drift on the power flux into the outer divertor is illustrated in figure 8. Solid lines show total, ion plus electron, convective plus conductive, power flux through the cell face boundary indicated by a thick red line in figure 6. This boundary is two rows of cells above another boundary separating 'divertor SOL' (using EDGE2D nomenclature) from the outer divertor. The power flux through it collects almost all radial power flux crossing the separatrix through the outboard (lower field) side, at the same time, the power crossing this boundary is not much affected by power sinks in the divertor and large cells surrounding the X-point. Dashed lines in figure 8 show only convective ion plus convective electron power fluxes through the same boundary. It is clear that convective power fluxes in the near SOL make a significant contribution to the total power fluxes. Convection is also responsible for the modulation of the shape of the total power flux profiles in the SOL. The modulation is clearly caused by the modulation of the shape of the E r profile shown in figure 2.
There is one extra factor supporting drift effects in the near SOL. In SOL rings closest to the separatrix, owing to the large flux expansion near the X-point field lines are longer, hence the parallel conductive power flux   potential distribution at the OMP is unaffected by long L || . Hence, the relative contribution of the convective power flux carried by the poloidal E×B drift is stronger than at other, more outward positions in the SOL. Another contribution towards the increase of a relative contribution of convection to total power fluxes is higher T i than T e , which is often the case in the 'main' SOL, upstream of the divertor. In the cases shown here the T i /T e ratio at the cell face boundary for which power fluxes are plotted in figure 8 is≈1.4, which increases the contribution of ion convective power fluxes.

Influence of transport coefficients on near SOL E r
As was pointed out in section 3, all cases were run with the same anomalous transport coefficients for all isotopes. Since particle and energy confinement, as well as transport coefficients at the plasma edge extracted from fluid 2D edge code simulations, indicate better confinement in plasmas with heavier hydrogen isotopes (see section 2), the question arises whether specifying the same transport coefficients for all isotope runs is realistic. The motivation for using the same transport coefficients was explained at the beginning of section 4. However, once the cause of isotope effects on the near SOL E r has been established, one can explore a possibility of the near SOL E r , via its effect on turbulence suppression, reducing transport coefficients. This effect could amplify isotope effects on the confinement. E r profiles for the isotope scan cases with variable transport coefficients are shown in figures 9(a)-(c). The D case is the original case, for which results were shown in previous figures. For H and T cases transport coefficients were changed according to the / m m D i scaling, where m i is either H or T mass. The square root dependence is chosen based on the results of EDGE2D-EIRENE simulations for JET L-mode plasmas, which established that in the D discharge D ⊥ in the SOL was by factor 1.4 lower in the D than in the H discharge [21]. Figure 9(a) shows the same profiles for original (same for all isotopes) transport coefficients, it only differs from figure 2 by the vertical scale which is adjusted to be the same as in figures 9(b), (c). Figures 9(b) shows results of varying only D ⊥ by factor / m m .  Figure 9(c) shows results of varying both D ⊥ and χ e,i by this factor. It is clear from figures 9(a)-(c) that reducing transport coefficients, mostly D ⊥ , leads to an increase in the near SOL E r . For the T cases, the E r spike, measured as the difference between the maximum E r value in the SOL and its value at the core ring closest to the separatrix, is larger in figure 9(c) than in figure 9(a) by factor 1.48. For the H cases, the same quantity is smaller by factor 1.42 in figure 9(c) than in figure 9(a). The near SOL E r spikes in these cases therefore appear to scale approximately as / m m , i D that is, inversely proportional to the variation in transport coefficients. The changes in E r are caused by changes in target T e profiles which in turn are caused by changes in transport coefficients.
Another test was made using only D cases. Figure 10 shows the effect of varying all transport coefficients by factors 2: starting with the original transport coefficients, they were multiplied by factors 2 and ½. A very strong impact of such a variation of transport coefficients on the near SOL E r spikes can be seen in the figure. Target T e profiles for the same cases as shown in figure 10, are shown in figure 11. Convective and total power fluxes near the entrance to the outer divertor for theses cases are shown in figure 12. In black (marked 'orig_case') are the same fluxes as shown in figure 8 for the D case, while in red (marked 'coef_×0.5') and blue (marked 'coef_×2') are D cases with transport coefficients multiplied by factors 0.5 and 2, respectively. The strongest impact of the variation of transport coefficients on E r profiles is seen in the first 4-5 rings in the SOL. Changes in E r values in the core are caused by changes in both electron density and ion temperature profiles: for the same outward particle and power fluxes n e and T i gradients, which determine E r in the core, vary inversely proportionally to transport coefficients.
Results of this sub-section show that positive feedback loops amplifying the influence of the near SOL E r on the E×B shear turbulence suppression may be considered, both for the formation of the E×B shear and for its effect on transport coefficients in the SOL. It has to be noted however that the modeling results presented here to a large extent depend on the existence of near SOL T e spikes, which are very likely to exist e.g. in L-mode plasmas close to the L-H transition. Where such spikes do not exist, for example, in detachment conditions, these arguments do not apply. Other conditions where the near SOL E r arguments should not apply, even though they must exist there, are probably H-modes with high input power, significantly exceeding the H-mode power threshold, where a very strong transport barrier inside of the separatrix across the pedestal region is already formed, with strong suppression of turbulence by the E×B shear.

Conclusions
Possible turbulence suppression by the E×B shear in the near SOL is likely to be an additional mechanism responsible for large scatter in experimental data on the H-mode power threshold (factor two in JET-ILW plasmas, depending on the divertor configuration), rather than be responsible for main trends in this threshold as described by power threshold scalings. Recent findings from AUG seem to have provided convincing evidence that scalings for the H-mode threshold power, in particular the Martin scaling [12], are related to the threshold in the ion power flux through the separatrix [28], which is not directly related to divertor target T e profiles and E r in the SOL. At the same time, the above findings do not explain the dependence of the H-mode threshold power (P LH ) on extreme edge region, as evidenced by the strong influence of the divertor configuration on P LH . It is therefore quite possible that the extreme edge, in particular SOL and divertor regions, can make a strong impact on P LH and account, at least partly, for a large scatter of experimental P LH (factor four scatter as seen in [12]).
Analysis of possible physical mechanisms for the influence of hydrogen isotope effects on the plasma particle and energy confinement in the plasma core is beyond the scope of Figure 10. Effect of varying all transport coefficients by factors 2 on E r at the OMP versus distance from the separatrix at the OMP position for the D case: starting with the original transport coefficients, they were multiplied by factors 2 and ½.  the present paper. It focuses on the extreme edge region, in particular, separatrix and near SOL regions, and a possible influence of the near SOL E r on the turbulence suppression in the SOL. EDGE2D-EIRENE modeling of hydrogen isotope scan cases shows that different neutral ionization patterns, depending on the isotope, may be responsible for the difference in the near SOL E r . It is well known that the largest neutral circulation takes place in the divertor, causing large ionization source there. The remaining neutrals, which were not ionized in the divertor, can penetrate (mostly) into the core and be ionized there. The role of CX neutrals coming from both the divertor and the main SOL in ionization sources deeper in the plasma must be very important. Large plasma volume around the X-point location presents a strong barrier for neutrals on their way to the core. The lightest isotope, hydrogen, has the largest probability to reach the X-point region and get ionized in the 1st poloidal ring, just outside of the separatrix. This ionization increases n e and reduces T e in this ring at the target, leading to lowest near SOL E r out of the three isotope cases.
The effect of the poloidal E×B drift may amplify the initial difference in the near SOL E r between EDGE2D-EIRENE cases for different isotopes caused by the difference in the ionization pattern. The extra convective poloidal power flux associated with the poloidal E×B drift causes higher target T e,i , which in turn result in larger near SOL E r . Such a self-amplifying effect may create larger difference in near SOL E r for cases with different isotopes compared to those caused only by the ionization patterns.
A potential positive feedback loop with a possible impact on the near SOL E r was identified in the modeling. In a series of cases in deuterium with different anomalous transport coefficients larger near SOL E r spikes are found in cases with lower transport coefficients. The larger near SOL E r in turn can result in turbulence suppression by the poloidal E×B shear. In relation to the isotope scan cases, variation of transport coefficients by factor / m 1 i led to the increase in the near SOL E r spikes by approximately factor m .
i Thus, again, a self-amplifying mechanism may be considered by which edge plasmas with heavier hydrogen isotopes may end up having stronger turbulence suppression and lower transport coefficients.