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

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²², 1.29 × 10²² and 1.06 × 10²² s⁻¹ 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.

Zoom In Zoom Out Reset image size

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⁻³ s⁻¹) 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.

Zoom In Zoom Out Reset image size

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²¹ s⁻¹, which has to be compared with the total ionization rate of 2.81 × 10²³ s⁻¹. 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 ${T}_{{\rm{e}}}^{7/2}/{L}_{| | }$ 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.

Zoom In Zoom Out Reset image size

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, $5/2{\rm{\Gamma }}T,$ where ${\rm{\Gamma }}$ 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: $\displaystyle \frac{{\rm{d}}T}{T{\rm{d}}t}=\displaystyle \frac{2}{3}\displaystyle \frac{{\rm{d}}n}{{\rm{d}}t}.$ 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.

Zoom In Zoom Out Reset image size

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 ${T}_{{\rm{e}}}^{7/2}/{L}_{| | }$ is lower. For example, for the 1st ring outside of the separatrix, the distance from the outer to inner target along the field line is 119.2 m, for the 3rd ring—109.2 m, for the 5th ring—96.6 m, and for the last, 20th ring—77.5 m. At the same time, electric 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.

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 $\sqrt{{m}_{D}/{m}_{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 $\sqrt{{m}_{D}/{m}_{i}}.$ 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 $\sqrt{{m}_{i}/{m}_{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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.