Numerical implications of including drifts in SOLPS-ITER simulations of EAST Open Access

Modeling of the plasma edge is widely used to understand which physical mechanisms can help reduce the power and particle exhaust toward the divertor targets in tokamaks.¹ For that, transport equations based on the ones of Braginskii are often solved. The SOLPS-ITER code package couples the finite volume plasma multi-fluid code B2.5 and the Monte Carlo kinetic neutral code EIRENE.^2,3 B2.5 solves the fluid transport equations for the plasma, and EIRENE calculates the source terms appearing in these equations resulting from interactions with neutrals.

Drift effects cause up-down pressure asymmetries and radial ion flows.⁴ This gives, in H-mode experiments, an asymmetry between the inner and outer divertor targets as observed experimentally in Ref. 5. In Refs. 6–9, it is shown for different tokamaks that drifts are required to obtain this asymmetry in SOLPS-ITER profiles: higher density and particle fluxes are observed at the inner target—the so-called high field side high-density region–, while higher temperatures and heat fluxes are observed at the outer target. To analyze further in detail the Ne-seeding experiments at EAST from which a first analysis was performed in Ref. 10, SOLPS-ITER simulations including drifts are required.

This paper focuses on the numerical implications of turning on drifts in SOLPS-ITER while physical interpretation of the results and comparison with experimental data will be dealt with elsewhere. As SOLPS-ITER is a 2D code, similar terms appear in the solved set of equations when drifts are introduced, making the code unstable. Therefore, several artificial numerical parameters are introduced in the set of equations solved in SOLPS-ITER.

The effects of these artificial numerical transport coefficients are verified. Additionally, the effects of space and time discretization and of the coupling to a Monte Carlo code are analyzed. It has been shown that these factors influence the outcome of SOLPS-ITER simulations from the numerical side and in that way might influence drifts-enabled simulations.^12–14 

The analysis is performed for the deuterium and Ne-seeded experiments from Ref. 10. Furthermore, the numerical implications of including Ne in an attached simulation are checked.

The performed analysis is limited to simulations for the EAST tokamak. For the first time, the effect of these numerical parameters is studied in detail for a SOLPS-ITER simulation of EAST. The findings are not necessarily valid for simulations of other devices, but as shown later on, some findings agree with published results elsewhere. Therefore, the presented outcome can help to obtain drift-enabled simulations for other devices.

The examined EAST discharges are in a disconnected double null configuration with the main upper divertor. However, as the separation between the two separatrices is $∼ 2 cm$⁠, they are considered as discharges in an upper single null configuration. If the upper single null topology is used during the grid generation for the SOLPS-ITER grids, this implies that the grids only extend $2 cm$ into the SOL outside the separatrix. As the power decay length for the examined experiments is estimated to be maximum $5 mm$⁠,¹⁰ a grid which goes $2 cm$ into the SOL covers 4 times the power decay length which should be sufficient.

The paper is organized as follows: in Sec. II, the way the electric currents are calculated in SOLPS-ITER is briefly discussed. In Sec. III, the effects of the choice of some numerical input parameters on the profiles of SOLPS-ITER simulations including drifts for attached purely deuterium, attached Ne-seeded simulations, and for the detached Ne-seeded case are presented. In Sec. IV, the impact of drifts on the distribution of the power deposition is examined, before coming to a summary of the main results.

In Ref. 11, the basic equations solved in the SOLPS-ITER code are given. They consist of the continuity equation, the parallel momentum conservation equation, and the ion and electron energy conservation equations. This set of equations can be solved while neglecting drift and electric current effects. To include drift and current effects on the particle flows in the simulations, in addition, the charge continuity equation has to be solved.

This shows that the plasma current consists of an (artificial) anomalous (⁠ $j A N$⁠) and a neoclassical current. The latter can be further split up in the diamagnetic current (⁠ $j ̃ dia$⁠), current caused by inertia and gyroviscosity (⁠ $j i n$⁠), viscosity contributions (⁠ $j ̃ vis | | + j ̃ visq + j vis ⊥$⁠), neutral-ion friction contributions (⁠ $j s$⁠), and a parallel current originating from the electron parallel momentum balance (⁠ $j | |$⁠).¹⁵ The tilde on some of the terms indicates that divergence free parts of these currents are eliminated analytically.¹⁶  $j ̃ dia , j ̃ vis | |$⁠, and $j ̃ visq$ represent physically a vertical guiding center drift caused by the curvature of the magnetic field B and $∇ B$⁠. Important here is to state the meaning of the artificial anomalous current. Often “anomalous” plasma parameters are referred to as a result of turbulence. In the presented paper, however, “(artificial) anomalous” is used to denote artificial transport quantities that are introduced to enhance the numerical stability and convergence of the code. This choice is made to use the same terminology as employed in the SOLPS-ITER manual.¹² 

The value of the neoclassical conductivity depends on the collision frequency regime in which the tokamak is operating. Therefore, it is often easier to ensure that $j A N ≪ j NEO$ at each location.

Similar to the artificial anomalous current, the artificial anomalous thermo-electric coefficient is only present for convergence purpose. Therefore, its contribution should be small in comparison with the classical thermo-electric coefficient: $α A N ≪ α C L$⁠.

To close the charge continuity equation [Eq. (1)], boundary conditions (BCs) should be imposed. At the far scrape-off layer, it is common to impose a zero gradient for $ϕ$⁠.^17,18 At the core boundary, the radial electric current is put equal to the diamagnetic current,^11,17 and at the private flux boundaries of the grid, the artificial anomalous current is put to zero. At the divertor targets, sheath BCs are imposed. As this results in a non-linear BC, a modified implementation of sheath BCs for drift simulations has been developed.¹² As mentioned in Ref. 5, the imposed sheath conditions limit the numerical stability of the code.

As the drift velocity will influence the overall velocity of the plasma particles, it also influences the other four transport equations (particle balance, parallel momentum balance, and energy equations for electrons and ions) as demonstrated in Ref. 11. Therefore, their computation has a large influence on the examined plasma quantities and is further investigated through this paper. For the studied EAST case, the typical directions of drift flows caused by the electric and magnetic field are indicated in Fig. 1.

Starting from the SOLPS-ITER setups of Ref. 19 for EAST, drifts are turned on. For the presented simulations, the same convergence criteria as in Ref. 19 are used. In the entire paper, the 3.0.7 master version of the SOLPS-ITER code is employed.

This is demonstrated in Fig. 2. The variations in plasma potential in the divertor region for a simulation in attached conditions are larger than the variations for a detached simulation. Therefore, it is expected that drifts are easier to include in a detached simulation than in an attached one.

Switching on all terms of the charge continuity equation at once does not lead to a converged solution. Whereas the viscous terms could be turned on fully straight away, the other drift terms (their strengths) have to be gradually added, step by step. First only $10 %$ of drift terms, afterward $50 % , 80 %$⁠, and finally $100 %$⁠. Only when the case with $10 %$ drifts is converged, the output is used to start the case with $50 %$ activated drifts, and so on.

The methods described in Ref. 14 are applied to the performed EAST simulations to enhance the inclusion of drifts: working with an intermediate solution by using a higher $σ ( A N )$ (and decreasing it again for the last simulation steps achieving final convergence), separating the time derivatives for the energy and density between inside and outside the separatrix, and including an artificial particle source to balance pumping and puffing more easily. Especially separating the time derivatives is necessary to obtain convergence for the presented simulations. More details about this method can be found in Ref. 14.

In this section, first the inclusion of drifts in attached purely deuterium simulations is discussed. Afterward, Ne-seeding is added, but the plasma is still attached. This changes some of the investigated numerics. Finally, enough Ne-seeding is included to bring the simulation in detachment. As indicated before, including drifts in detached cases is easier than in attached ones.

It is verified how the process of including drifts can be done in the most efficient way for a purely deuterium SOLPS-ITER simulation of EAST. The investigated simulation is the reference case from Refs. 10 and 19. Equations (7) and (8) show that the computation of drift velocities is influenced by two factors: the $ϕ$ and B fields at one hand and the calculation of their derivatives on the other hand. First, the influence of the imposed artificial anomalous conductivity and thermo-electric coefficient on the $ϕ$ field is investigated. As these artificial anomalous terms affect the solved set of equations, they modify the investigated physical problem. Afterward, the influence of the discretization is investigated by the effect of the grid on the computation of $ϕ$⁠. In a next step, the influence of the number of EIRENE particles and the averaging over simulation steps is discussed, before finishing with an analysis of the used time step and numerical switches in SOLPS-ITER. For all presented purely deuterium simulations, an input power of $2.05 MW$ is imposed at the core boundary of the grid.

Instabilities in the solution of the potential equation might lead to artificial oscillations or spikes in the resulting potential field which is determined by the set of solved equations. Therefore, the influence of the choice of the artificial anomalous parameters—σ_AN and α_AN—on the potential is investigated.

The artificial anomalous parameters, which were successfully used to include drifts in ASDEX Upgrade,¹⁷ are used as a starting point and adjusted to the needs of the presented work. This results in $c σ α , 0 = 5.0 × 10 − 5$ and $c α , 0 = 1.0 × 10 − 5$⁠. In Fig. 3, it is verified that $σ A N ≪ σ NEO$ by examining that $j ̃ A N ≪ j ̃ NEO$ for the simulations discussed above at the outer midplane. Similar results are obtained at the upper outer target (UOT). This confirms that $c σ α , 0 = 5.0 × 10 − 5$ and $c α , 0 = 1.0 × 10 − 5$ are appropriate choices for the performed EAST simulations.

First, the impact of σ_AN on $ϕ$ is investigated; afterward, the impact of α_AN is studied.

While changing the values of $c σ α , 0$⁠, the impact on the convergence becomes clear. Where a simulation in which $c σ α , 0$ is increased to $10.0 × 10 − 5$ converges within a couple of days, it was not possible to get a converged solution in case of $c σ α , 0 ≤ 4.0 × 10 − 5$ as the solution for the density and temperature profiles kept oscillating and did not become stable.

In Fig. 4, the plasma potential is compared at the upper outer target (UOT) for different values of $c σ α , 0$⁠. While the convergence speed increases with increasing $σ α , 0$⁠, artificial smoothing of the potential target profile is observed. This indicates a smaller gradient of the target potential, and in that way, a smaller electric field. This means that we have to make a trade-off between accuracy and convergence speed.

The potential drop of Fig. 4 causes opposite poloidal E × B drift velocities across the separatrix at the UOT as shown in Fig. 5.

As the ion saturation current j_s can be compared with the divertor Langmuir probes at EAST, this is one of the key plasma quantities for comparison with experiment.¹⁰ The effect of the employed artificial anomalous conductivity on the ion saturation current at the UOT is shown in Fig. 6: for the simulation with $c σ α , 0 = 10.0 × 10 − 5$⁠, the peak value for the ion saturation current is $∼ 50 %$ lower than for the simulation with $c σ α , 0 = 4.5 × 10 − 5$⁠, but Fig. 4 shows that this is because of a smaller gradient in the potential.

In Fig. 6, the simulation result without drift effects is displayed in purple to show the influence of drifts on the j_s profile. In red, the measurements from the divertor Langmuir probes for the purely deuterium plasma discussed in Ref. 10 are plotted. This indicates that the inclusion of drifts in the simulation improves the agreement between simulations and experiments which is the motivation for the effort presented in this paper to include drifts in SOLPS-ITER simulations. The figure indicates that the peak value for j_s decreases for larger $c σ α , 0$⁠. To get an overall quantification of the improved agreement between experimental and simulation data, the norm of the difference between them is calculated. Using the larger value for $c σ α , 0$ decreases this norm with $16.9 %$⁠. If the peak value for j_s is not included in the calculation of the norm, only an improvement of $1.65 %$ is obtained. This indicates that mainly one value is changing due to the increased $c σ α , 0$⁠. As indicated before, the observed changes originate from the artificial smoothing of the potential target profile and do not improve the described physics.

As the grid cell width increases from the outer midplane (OMP) toward the targets due to the flux expansion, this means that the artificial anomalous conductivity is increased in the region around the targets. In this region, the effect will be similar to increasing $c σ α , 0$ as demonstrated before, but in the remaining part of the SOL, the effect will be smaller. Figure 7 shows, indeed, that the potential drop is eliminated using $c σ α , 7 = 1$⁠. Also in this case, it is easier to achieve convergence, but there is a loss of accuracy. The effect on the j_s profile at the UOT is shown in Fig. 8. It might be suggested this method is better than increasing σ_AN everywhere as only the divertor region is affected.

In a next step, the influence of α_AN is investigated. In Ref. 12, it is stated that at least σ_AN or α_AN cannot be zero if drift flows are included in a SOLPS-ITER calculation. As σ_AN is fixed now to a finite value, the effects of putting $c α , 0 = 0.0 × 10 − 5 , c α , 0 = 1.0 × 10 − 5$⁠, and $c α , 0 = 2.0 × 10 − 5$ are examined. The resulting potentials are shown in Fig. 9. Also here, the simulations will converge faster for higher values of $c α , 0$⁠, but the gradient in the potential is smaller and, thus, the peak in j_s is smaller. The effect on the ion saturation current at the UOT is shown in Fig. 10.

It should be remarked that the investigated artificial anomalous quantities have the largest influence on the E × B drifts. The diamagnetic drift velocities mainly depend on the magnetic field which is an input parameter for SOLPS-ITER and independent of the potential. Therefore, their values stay similar as indicated in Fig. 11 for the diamagnetic drift profiles at the UOT for different values of $c σ α , 0$⁠.

In Sec. III A 1, it was shown that the artificial anomalous conductivity and thermo-electric coefficient can speed up convergence but decrease the accuracy of the calculated solution.

In this section, the effect of the grid on the calculated solution is investigated.

Figure 12 shows the potential at the UOT and Fig. 13 the ion saturation current at the UOT in case of the original grid of Ref. 19 (original grid) and in case of a grid with the same number of grid cells, but with a poloidal and radial grid width of minimum 1 mm (modified grid). Normally, SOLPS-ITER grids are refined toward the divertor targets as the plasma–wall interaction has a major influence on the simulation result and as you need to resolve larger gradients with respect to the rest of the simulation domain. In the modified grid, this refinement is not performed as can be seen in Fig. 14. In that way, grid cell widths smaller than 1 mm are avoided. The modified grid is the one used in all presented studies in this paper. As can be seen from Fig. 14, both grids cover the same region in EAST.

While there are numerical oscillations in the potential calculation using the original grid, they are not present using the modified one. In the original grid, the minimal radial cell width is 1 mm and the poloidal one minimum 0.667 mm. In Fig. 15, the energy crossing the separatrix is investigated for a further refined grid of Ref. 19 (refined grid) in which also the radial cell width is below 1 mm. The figure shows that for such a refined grid, a numerical instability occurs in the simulation. This is in agreement with the findings of Ref. 18 where a similar instability appeared at the divertor targets for TCV simulations with radial grid cell widths below 1 mm.

Increasing further the minimal poloidal and radial grid width results in even smoother potential profiles. However, as this leads to a larger discretization error, this is not investigated in detail. Coarsening the grid increases in that way the numerical diffusion which gives a similar effect as an increase in $c σ α , 0$ and $c α , 0$ discussed in Sec. III A 1. In the results presented in this paper, the modified grid of Fig. 14(b) with a radial and poloidal grid cell width of minimum 1 mm is employed.

Next, the influence of the number of involved Monte Carlo particles in the EIRENE simulations is investigated. The default number of EIRENE particles included in the presented simulations is based on Table 1 of Ref. 19, and the parameters used in the presented study are reported in Table I. After the initial transient, a simulation is averaged over 45 000 iterations to decrease the statistical error. In Ref. 19, it is shown that this combination of the involved number of Monte Carlo particles and SOLPS-ITER averaging results in statistical errors below $0.35 %$ on all investigated plasma quantities for a simulation without drifts. In Fig. 16, the effect of including 10 times more or 10 times less Monte Carlo particles in the simulation with fully activated drifts is demonstrated. This indicates that the reduced fluctuations between different iterations—resulting from an increased number of Monte Carlo particles—influence the potential profile.

Figure 16 supports the findings of Sec. III A 2, namely, that a higher potential drop around the separatrix is more realistic. It also indicates that with the default number of EIRENE particles, the potential drop around the separatrix is still underresolved. It should be noted that the analysis of Ref. 19 showed that for a simulation without drifts, this number of EIRENE particles is sufficient.

Figure 17 shows the importance of SOLPS-ITER averaging over several simulations.²³ The figure shows that a larger potential drop is found for the non averaged simulations, but this is caused by a deeper drop around the separatrix, while Fig. 16 and the analysis of Sec. III A 2 show that this is not expected (the potential should go higher before the separatrix, not lower afterward). It should be remarked that the profiles can vary from simulation to simulation if no averaging is applied. On top, also the profiles at other locations are not as smooth as in the averaged case.

The next topic discussed is the required time step for the simulations. Where Secs. III A 1, III A 2, and III A 3 showed influences on the final simulation result, the time step rather impacts the convergence itself. For the simulations without drifts from Ref. 19, a time step of $5.0 × 10 − 5 s$ is required to obtain a stable solution. However, this is not possible for simulations in which the viscous terms are activated. In Fig. 18, the variations in time of the total energy crossing the separatrix (a) and core boundary (b) are given for simulations in which the viscous terms are activated. As can be seen in yellow on Fig. 18(b), with a time step of $1.0 × 10 − 6 s$⁠, the simulation does not converge and gives numerical oscillations for the energy crossing the core boundary. A time step of $1.0 × 10 − 7 s$ is needed (the blue line in both plots of Fig. 18). The requirement of a small time step can be explained by the segregation of all different equations and the intrinsic noisy interaction with the EIRENE Monte Carlo code. This requires small time steps to avoid numerical instabilities in drift cases as the maximum time step is related to the convergence of the process with the smallest timescale involved in the numerical problem.¹⁸ In SOLPS-ITER, this means that a higher time step leads to divergence in the solution of the potential equation.¹⁴ 

In this process, increasing the value of $b c ref , t e / t i$—called boundary under-relaxation factors in the further scope of this paper—means less under-relaxation and is the same as putting a very strong forcing toward the desired flux $Q imposed , e / i$⁠. $Q e , old$⁠, and $Q i , old$ are the actual power fluxes going through the core boundary in the previous time step, where $Q imposed , e$ and $Q imposed , i$ are the ones imposed in the boundary condition. This effect is demonstrated with the red curve in Fig. 18: if you try to force the simulation to reach the wanted power flux too fast, the simulation will not converge. The fact that the influence of the previous time step must be limited to avoid divergence also contributes to the long converging times required for drift simulations.

Additionally to the imposed under-relaxation factors, also an appropriate choice of other numerical parameters is needed. In a simulation in which no drifts are present, the plasma equations are solved with a time step 100 times higher inside the separatrix for the ion and electron energy equations, and a time step 10 times higher for the continuity equation. As the critical point in drift simulations lies inside the separatrix, this speed-up method does not lead to a converged simulation anymore. Therefore, the time step inside the separatrix should be smaller than the one outside the separatrix as discussed in Ref. 14. In the presented work, the time step inside the separatrix is taken 100 times smaller than the one outside the separatrix.

In a next step, Ne-seeding is added to SOLPS-ITER simulations. As neon will also radiate in the core, the input power crossing the core boundary of the SOLPS-ITER grid is decreased from $2.05$ to $1.80 MW$ based on the bolometer data for the investigated discharge.¹⁰ Other boundary conditions are kept the same as in the unseeded simulation. If only a small amount of Ne is added, the simulations are still attached. This makes it possible to examine the influence of the presence of impurities on the numerics in the SOLPS-ITER code. More precisely, the time step and numerical parameters influencing the convergence are examined in the present work.

In Fig. 18, it was demonstrated that a decreased boundary under-relaxation factor is required to converge the simulations. Figure 19 shows the opposite effect if Ne-seeding is present: a too low under-relaxation factor decreases the influence of the current time step too much. This makes that the simulation does not converge to the correct energy crossing the core boundary of the employed grid. In case the under-relaxation factor is increased, the power crossing the core boundary increases up to $98 %$ of the power which is imposed at the core boundary. If the boundary under-relaxation factor is increased too much, the simulation converges faster, but there is a discrepancy between the imposed power and the result of the simulation. Figure 19 shows that a boundary under-relaxation factor of 0.03 (increased under-relaxation) still gives a correct power crossing the core boundary. In Fig. 20, the effect of the boundary under-relaxation factor on the instant ion saturation current at the UOT is shown for the Ne-seeded simulation. This indicates that a larger boundary under-relaxation factor leads to lower peaks as the set of equations is numerically more stable. Therefore, a boundary under-relaxation factor of 0.03 (instead of 0.001 in the unseeded, purely deuterium simulations) is used for the performed Ne-seeded simulations.

It should be remarked that these simulations were performed with 10 times less EIRENE particles than the unseeded simulations presented earlier in this paper. Figure 16 shows that the potential in an unseeded simulation is smoother if less EIRENE particles are used.

Figure 2 showed a decreased potential in detached simulations. Nevertheless, the potential is not completely smooth as demonstrated in Fig. 21 with the potential profile at the UOT in a simulation in which enough Ne is added to bring the plasma in detachment. Also for this simulation, 10 times less EIRENE particles are used due to the limited available computing time. In contrast to the earlier shown simulations, temperature boundary conditions are applied at the innermost grid boundary. This means that no under-relaxation factors are needed. However, this makes it more difficult to mimic the physics of a real experiment as the input power is known more precisely than the temperature at the core boundary of the B2.5 grid.

As the grids in the 3.0.7 master version of SOLPS-ITER are not extended up to the first wall, it should be verified that they cover the region of the tokamak edge where the most significant fraction of physical processes take place. Reference 24 showed for JET that, if the grid of the simulation is too narrow, not all the power is captured in the simulation. In order to determine if the grid for the examined EAST discharges extends far enough in the SOL to capture most of the power in the simulation, it is important to investigate how the injected power is dissipated. If too much power is deposited on the B2.5 grid boundary toward the main chamber wall, the simulation will not be able to give an accurate description of the power deposition in the divertor. Therefore, the power balance is verified. Special attention is put on the effect of drifts on the power crossing the grid boundary closest to the first wall as indicated in Fig. 22. For the examined simulations (with and without drifts), the power crossing this grid boundary is summarized in Table II. As it was found that the number of EIRENE particles does not influence the power balance, all percentages summarized in Table II are given for simulations with 10 times less particles as in that case simulations for all examined conditions were available.

This overview demonstrates that in all simulations $∼ 75 %$ of the energy crossing the core boundary of the grid (in red in Fig. 22) is not deposited on the first wall. In the deuterium simulations, most of the energy will be transported toward the inner and outer divertor legs. In the Ne-seeded simulations, on the other hand, a significant part of the power will already be radiated before the entry of the divertor. The power sharing between inner and outer divertor shown in Table II clearly indicates that drifts are required to obtain an asymmetry between inner and outer divertor. In the presence of drifts, only a small fraction of the power is going toward the inner divertor, and most of the power enters the outer divertor leg. This is also nicely illustrated on the 2D plots of Fig. 23, which show the 2D profiles for the density and temperature with and without drifts. In case drifts are included, a clear asymmetry between the inner and outer divertor is present while this is absent when drifts are not switched on. As indicated before, this is in agreement with findings on other devices.^6,7

The remaining $∼ 25 %$ of the energy transported from the core into the scrape-off layer will be dissipated volumetrically outside the B2.5 grid and eventually mainly deposited on the first wall and is not considered further in the simulation. A simulation with a grid extended until the first wall as introduced in Ref. 25 for SOLPS-ITER would make it possible to investigate where this energy is exactly deposited. However, as the simulation domain for the current simulations covers $73.9 %$– $80.5 %$ of the power deposition, they can all (both without and with drifts) be used to investigate more in detail the power deposition in the divertor of EAST.

In this paper, the numerical implications of including drifts in the SOLPS-ITER simulations of EAST are studied. For that, three simulations in upper single null configuration are analyzed: an attached purely deuterium simulation, an attached simulation in which low levels of Ne-seeding are present, and a detached Ne-seeded simulation. As the gradients in potential profiles strongly influence the plasma flows, and as these potential gradients are larger during attachment, the attached drift simulations are more difficult to converge than the detached ones. Therefore, also the influence of employed numerical parameters in the SOLPS-ITER setup is larger.

In order to converge SOLPS-ITER drift simulations, artificial anomalous terms are introduced in the charge continuity equation. These terms are determined by an artificial anomalous conductivity and an artificial anomalous thermo-electric coefficient. In the presented work, it is shown that the choice of these parameters largely influences the potential, and in that way, the $E × B$ plasma flows. The diamagnetic drifts, on the other hand, are independent of these artificial anomalous parameters.

In III A 2, the influence of the grid on the plasma potential is verified. It came out that grid cell widths below 1 mm give rise to numerical instabilities which are reflected in the potential profile. Therefore, it is recommended to use grid cell widths of at least 1 mm in the radial and poloidal direction.

The number of EIRENE particles in the SOLPS-ITER simulation also influences the potential profile. Around the separatrix, an abrupt change in the plasma potential is expected. However, this peak can give numerical issues as the Monte Carlo noise of the EIRENE simulation resulting from too few particles gives rise to local fluctuations which may impact convergence. Therefore, a larger number of EIRENE particles is necessary. This has as the drawback that the required computing time increases. Additionally, applying SOLPS-ITER averaging after convergence is reached and smooths the potential profiles.

Next, the time step and involved under-relaxation factors for the imposed energy boundary conditions which are needed to obtain a converged simulation are described.

After finishing the analysis for the attached purely deuterium simulation, the attached Ne-seeded and detached Ne-seeded simulations are analyzed. The choice of under-relaxation factors for the energy boundary condition seems to be dependent on the specific case: the attached Ne-seeded simulation only converged correctly if increased under-relaxation factors were imposed.

In the last section, the effect of drifts on the power crossing the grid boundary toward the main wall is investigated. As the employed SOLPS-ITER grid is not extended up to the first wall, it is important to keep track of this to avoid an artificial loss of power through the outermost boundary of the B2.5 grid. The performed analysis shows that the simulations can be used for divertor power deposition studies of the investigated EAST discharges.

The authors would like to thank the EAST team, in particular Yunfeng Liang, Liang Wang, Fang Ding, and Kedong Li, for providing the EAST-related parameters (equilibrium, geometry, pump efficiency…) used in the present study.

This work has partially been carried out within the EU-CN collaboration framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

This work was funded in parts by the U.S. Department of Energy (No. DE-SC0014210) and the College of Engineering at UW Madison, WI, USA.

The SOLPS-ITER simulations were performed at the Marconi supercomputer from the National Supercomputing Consortium CINECA.

The authors have no conflicts to disclose.

D. Boeyaert: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Methodology (equal); Software (equal); Validation (lead); Visualization (lead); Writing – original draft (lead). S. Carli: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Supervision (supporting); Validation (equal); Writing – review & editing (lead). W. Dekeyser: Methodology (equal); Software (equal); Supervision (supporting); Writing – review & editing (equal). S. Wiesen: Funding acquisition (equal); Project administration (equal); Supervision (equal). M. Baelmans: Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.