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

The inclusion of drifts in plasma edge codes like SOLPS-ITER is required to match simulation data with experimental profiles. However, this remains numerically challenging. In this paper, the effect of some numerical factors on the final plasma solution is investigated. This study is performed on three EAST simulations in upper single null configuration: an attached purely deuterium case, an attached case with limited Ne-seeding and a detached Ne-seeded case. The effects of the anomalous conductivity and anomalous thermo-electric coefficient on the plasma potential are investigated. Next, the effect of the employed grids is shown. In order to investigate these effects, accurate drift simulations are needed. Therefore, the employed time step and used numerical parameters are discussed for the three studied simulations. For all presented simulations, it is verified that the restriction of the grid to the first flux surface tangent to the main chamber wall has only a small effect on the divertor solution. This means that the main power dissipation takes place inside the simulated domain and only a small fraction of the power is leaving the B2.5 grid trough the grid boundary closest to the first wall. Finally the effect of drifts on the asymmetry between the inner and outer target for EAST simulations is demonstrated.


I. INTRODUCTION
Modeling of the plasma edge is widely used to understand which physical mechanisms can help reducing the power and particle exhaust towards the divertor targets in tokamaks 1 .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 4 .This gives, in H-mode experiments, an asymmetry between the inner and outer divertor targets as observed experimentally in ref. 5 .7][8][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 indicated in ref. 11 the numerical schemes in SOLPS-ITER become unstable for drift simulations when an H-mode transport barrier is imposed in the simulation.Therefore, several artificial numerical parameters are introduced in the set of equations solved in SOLPS-ITER.
The effects of these artificial numerical transport coefficients is 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][13][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 limits 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.The findings are not necessary 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 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 10 , 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 section 2, the way the electric currents are calculated in SOLPS-ITER is briefly discussed.In section 3, the effects of the choice of some numerical input parameters on the profiles of SOLPS-ITER simula-tions including drifts for attached purely deuterium, attached Ne-seeded simulations, and for the detached Ne-seeded case are presented.In section 4, the impact of drifts on the distribution of the power deposition is examined, before coming to a summary of the main results.

II. MODEL DESCRIPTION
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 effects on the particle flows in the simulations, the charge continuity equation has to be additionally solved.
The charge continuity equation (∇ • j = 0) is formulated as follows in SOLPS-ITER (using the toroidal symmetry of tokamaks): in which x is the poloidal direction and y the radial direction normal to the magnetic flux surfaces.h x is the cell width in the poloidal direction, h y the cell width in the radial direction.g indicates the surface area between a grid cell and its neighbor.In this equation j is the plasma current vector and it is split in the following way: This shows that the plasma current consists of an (artificial) anomalous ( jAN ), and a neoclassical current.The later can be further split up in the diamagnetic current ( jdia ), current caused by inertia and gyroviscosity (j in ), viscosity contributions ( jvis|| + jvisq + j vis⊥ ), neutral-ion friction contributions (j s ), and a parallel current originating from the electron parallel momentum balance (j || ) 15 .The tilde on some of the terms indicates that divergence free parts of these currents are eliminated analytically 16 .jdia , jvis|| , and jvisq 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 12 .
In modeled H-mode discharges with plasma edge codes like SOLPS-ITER, typically a transport barrier is imposed for the perpendicular diffusion coefficient 17 .This results in small values for the imposed perpendicular diffusion for a large fraction of the simulation domain.As a result, solving equation 1 is difficult from a numerical viewpoint.Therefore, the artificial anomalous current jAN is introduced which is dependent on the gradient of the potential.
This artificial anomalous current is not physical, but is only present to ensure smooth convergence.As a result, its contribution should be large enough to ensure convergence but at the same time small to avoid a non-physical simulation result.Therefore, an appropriate choice of the artificial anomalous conductivity σ AN should be made.This means that σ AN << σ NEO ≈ µ i1 (BR) 2 with µ i1 the parallel viscosity 14 .In SOLPS-ITER, σ AN is determined by c σ α,0 with the following formula 12 : in which n e is the electron density at the innermost grid surface at the outer midplane.Further possible dependencies of σ AN are described in ref. 12 but are not considered in the scope of the presented work.
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 AN << j NEO at each location.
As indicated earlier, the last right-hand term of equation 2 (j || ) is originating from the electron parallel momentum balance.If no impurities are present in the plasma, this current is given by the following equation 12 : in which σ || is the Spitzer conductivity and φ the plasma potential.The first part of the right hand side compensates for the diamagnetic current, and forms the Pfirsch-Schlüter current.The second term on the right is the thermo-electric current which arises due to temperature gradients.In this expression, α ex is the thermo-electric coefficient.This coefficient consists of a classical part, and an artificial anomalous contribution.As for the artificial anomalous current, the artificial anomalous thermo-electric coefficient is: Similar to the artificial anomalous current, the artificial anomalous thermo-electric coefficient is only present for convergence purposes.Therefore, its contribution should be small in comparison with the classical thermo-electric coefficient: α AN << α CL .
To close the charge continuity equation (equation 1), boundary conditions (BCs) should be imposed.At the far scrape-off layer, it is common to impose a zero gradient for φ 18,19 .At the core boundary, the radial electric current is put equal to the diamagnetic current 11,18 , 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 12 .As mentioned in ref. 5 , the imposed sheath conditions limit the numerical stability of the code.
Once the plasma potential and the magnetic field are known, the drift velocities resulting from the E x B drift and the diamagnetic drift velocities resulting from the ∇B drift can be calculated: 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 directions of the appearing drift flows caused by the electric and magnetic field are indicated in figure 1.
Starting from the SOLPS-ITER setups of ref. 20 for EAST, drifts are turned on.For the presented simulations the same convergence criteria as in ref. 20 are used.In the entire paper the 3.0.7 master version of the SOLPS-ITER code is employed.

III. INCLUDING DRIFTS IN SOLPS-ITER EAST SIMULATIONS
Equation 5 can be rewritten as: in which . Under attached conditions, the Spitzer conductivity will be high 21 meaning that the first right-hand term of equation 9 is small.The basic two-point model in attachment shows that 2n t T t = n u T u in which 'u' refers to the upstream conditions and 't' to the target conditions 22 .This makes that also the last right-hand term is small, resulting in: This shows that the local potential is determined by the integration of the temperature along a flux tube, and that the potential at a specific location is directly linked to the electron temperature at that location.Under attached conditions, the maximum temperature at a divertor target is larger than in detachment and the temperature gradients in the perpendicular direction are significant.As a consequence, large perpendicular gradients in the plasma potential will be present.This The direction of the (ion) drift flows in the studied EAST discharge.Three different flows are indicated: the one caused by the magnetic field and its derivative (in black), the one caused by the radial electric field in combination with the magnetic field (in red) and the one caused by the poloidal electric field and the magnetic field (yellow).The influence of the radial electric field is larger than the one of the poloidal electric field.
makes the charge continuity equation more unstable from a numerical viewpoint.This is demonstrated in figure 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, afterwards 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 σ (AN) (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.Afterwards 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.

A. Including drifts in a purely deuterium attached simulation
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,20.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.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.Afterwards, the influence of the involved numerics is investigated by the effect of the grid on the computation of φ .In a next step, the influence of the involved 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.

Appropriate choice of artificial anomalous conductivity and thermo-electric coefficient
A smooth profile of the plasma potential φ is required as discontinuities can hamper gradient calculations in equation 7. φ 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 18 , are employed in the presented work.This results in c σ α,0 = 5.0 • 10 −5 and c α,0 = 1.0 • 10 −5 .This choice is explained by the fact that ASDEX Upgrade has a similar size as EAST.In figure 3 it is verified that σ AN << σ NEO by examining that jAN << jNEO for the simulations discussed above at the OMP.Similar results are obtained at the UOT.First the sensitivity of σ AN on φ is investigated, afterwards 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 in which the solution for the density and temperature profiles kept oscillating and did not became stable.
In figure 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 an underestimation of the electric field.This means that we have to make a trade-off between accuracy and convergence speed.As the ion saturation current can be compared with the divertor Langmuir probes at EAST, this is one of the key plasma quantities for comparison with experiment 10 .The effect of the employed artificial anomalous conductivity on the ion saturation current at the UOT is shown in figure 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 figure 4 shows that this is because of an underestimation of the potential.
In figure 6 the simulation result without drift effects is shown in purple to show the influence of drifts on the j s profile.In red the measurements from the divertor Langmuir probes 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.
In section III A 2 the effect of small grid cells on the stability of the solution is demonstrated.A way to decrease the effect of such instabilities appearing with small grid cells, is changing σ AN depending on the grid cell width.In the 3.0.7 master version of SOLPS-ITER, the following scaling of σ AN with the grid cell width is available if c σ α,7 ̸ = 0: As the grid cell width is smaller around the outer midplane (OMP) due to the flux expansion towards the targets, 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 figure 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 figure 9. Also here, the simulations will converge faster for higher values of c α,0 but the potential profile is less well representing the drop around the separatrix.The effect on the ion saturation current at the UOT is shown in figure 10.It should be remarked that the investigated artificial anomalous quantities have the largest influence on the E x 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 figure 11 for the diamagnetic drift profiles at the UOT for different values of c σ α,0 .

of the spatial discretization
In the previous section it was shown that the artificial anomalous conductivity and thermo-electric coefficient can In this section, the effect of the grid on the calculated solution is investigated.
Figure 12 shows the potential at the UOT and figure 13 the ion saturation current at the UOT in case of the default grid of ref. 20 (default 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 towards the divertor targets as the plasmawall interaction has a major influence on the simulation result.In the modified grid, this refinement is not performed as can be seen in figure 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 figure 14, both grids cover the same region in EAST.While there are numerical oscillations in the potential calculation using the default grid, they are not present using the modified one.In the default grid, the minimal radial cell width is 1 mm and the poloidal one minimum 0.667 mm.In figure 15 the energy crossing the separatrix is investigated for a further refined grid of ref. 20 (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. 19 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 facilitates 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 c σ α,0 and c α,0 discussed in the previous section.In the results presented in this paper, the modified grid with a radial and poloidal grid cell width of minimum 1 mm is employed.

Influence of the number of included Monte Carlo particles and averaging
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. 20 and the parameters used in the presented study are copied in table I.After convergence, a simulation is averaged over 45,000 iterations to decrease the statistical error.In ref. 20 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 figure 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 particlesinfluence the potential profile.FIG.16: The influence of the number of involved EIRENE particles on the plasma potential at the UOT.
Figure 16 supports the findings of the previous sections, 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 underestimated.It should be noted that the analysis of ref. 20 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 23 .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 figure 16 and the analysis of the previous section show that this is not expected (the potential should go higher before the separatrix, not lower afterwards).On top, also the profiles at other locations are not as smooth as in the averaged case.

Employed time step and used numerical parameters to obtain convergence
The next topic discussed is the required time step for the simulations.Where the previous sections showed influences on the final simulation result, the time step rather impacts the convergence itself.For the simulations without drifts from ref. 20 , a time step of maximum 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 figure 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 figure 18b, 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 grid boundary.A time step of maximum 1.0 • 10 −7 s is needed (the blue line in both plots of figure 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 time scale involved in the numerical problem 19 .In SOLPS-ITER this means that a higher time step leads to divergence in the solution of the potential equation 14 .
In ref. 14 the appearing numerical instabilities are explained by "the poloidal redistribution of particles inside the separa- FIG. 18: The instabilities appearing when the viscous terms are included.In (a) the instabilities on the energy crossing the separatrix for different setups are shown and in (b) the ones while crossing the core boundary of the grid.'Original' in the legend refers to the default values for the boundary under-relaxation factors in SOLPS-ITER (bc re f ,te = bc re f ,ti = 0.01).For the blue curves, these values are lowered to bc re f ,te = bc re f ,ti = 0.001.
trix by E x B drifts in combination with the modification of the radial electric field by diamagnetic currents".These observations are obtained for ASDEX Upgrade and ITER simulations.As indicated in figure 18, a similar instability happens in the studied EAST discharges.The energy balance is examined during all simulations indicating that the energy crossing the core boundary of the grid does not stay stable.For the performed drift simulations, the energy crossing the core boundary is imposed as a boundary condition for the ion and electron energy equations.As such a flux BC is from the numerical viewpoint less stable than a Dirichlet type BC, drift simulations in SOLPS-ITER require a type of BC in which the code converts these imposed heat fluxes to electron and ion temperatures 12 .During a series of iterations, the new imposed temperature will depend on the actual fluxes imposed at the BC, the ones at the previous time step, and on how the imposed temperature should be corrected to match the heat flux from the BC: In this process, the importance of the previous time step can be regulated with under-relaxation factors -bc re f ,te and bc re f ,ti , called boundary under-relaxation factors in the further scope of this paper.f h e and f h i are the actual power fluxes going through the core boundary, where f imposed,e and f imposed,i are the ones imposed in the boundary condition.This effect is demonstrated with the red curve in figure 18: if the influence of the previous time step is too large, the simulation will not converge.The fact that the influence of the previous time step is limited, also contributes to the long converging times required for drift simulations.
Additionally to the imposed under-relaxation factors, also an appropriate choice of the 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 lower than the one outside the separatrix as discussed in ref. 14 .

B. Including drifts in a Ne-seeded attached simulation
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 MW to 1.80 MW.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 presented work.
In figure 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, however, 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 still gives a correct power crossing the core boundary.In figure 20 the effect of the boundary under-relaxation factor on the 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 is used for the performed Neseeded 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   Figure 2 showed a decreased potential in detached simulations.Nevertheless, the potential is not completely smooth as demonstrated in figure 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.

IV. EFFECT OF DRIFTS ON THE POWER CROSSING THE GRID BOUNDARY TOWARDS THE MAIN WALL
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.Ref. 24 showed for JET that, if the grid of the simulation is too narrow, it is difficult to draw physical conclusions as the decay-length type BCs impact the solution too much and the plasma state is determined by these BCs.In order to determine if the grid extends far enough in the SOL, it is important to investigate how the injected power is dissipated.If too much power is deposited on the B2.5 grid boundary towards 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 figure 22.For the examined simulations (with and without drifts) the power crossing this grid boundary is summarized in table II.As it came out 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 figure 22), is not deposited on the first wall.In the deuterium simulations, most of the energy will be transported towards 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 towards the inner divertor, and most of the power enters the outer divertor leg.This is also nicely illustrated on the 2D plots of figure 23, which show the 2D profiles for the density and temperature with and without drift.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 % to 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.V. SUMMARY 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 are 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 influence the potential, and in that way the E x B plasma flows.The diamagnetic drifts, on the other hand, are independent of these artificial anomalous parameters.
In a next section, the influence of the grid on the plasma potential is verified.It came out that grid cell widths below 1 mm give raise 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 a drawback that the required computing time increases.Additionally, applying SOLPS-ITER averaging after convergence is reached, 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 Neseeded simulations are analyzed.The choice of underrelaxation factors for the energy boundary condition seems to be dependent on the sepcific case: the attached Neseeded simulation only converged correctly if increased under-relaxation factors were imposed.
In a last section, the effect of drifts on the power crossing the grid boundary towards 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 ar-tificial 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.

FIG. 2 :
FIG. 2: The plasma potential for a simulation in the high-recycling regime (a) and the one for a detached simulation (b)

FIG. 4 :
FIG. 4:The influence of the choice of c σ α,0 on plasma potential at the UOT.

FIG. 5 :
FIG. 5: The influence of the choice of c σ α,0 on the V ExB x drift velocity at the UOT.

FIG. 6 :
FIG.6:The influence of the choice of c σ α,0 on the ion saturation current at the UOT.For reference the simulation result without drift effects is shown in purple, and the experimental data from the divertor Langmuir probe are given in red.

FIG. 7 :
FIG. 7:The influence of the choice of c σ α,7 on plasma potential at the UOT.

FIG. 8 :
FIG.8:The influence of the choice of c σ α,7 on the ion saturation current at the UOT.

FIG. 9 :FIG. 10 :
FIG. 9:The influence of the choice of c α,0 on plasma potential at the UOT.

FIG. 11 :
FIG.11:The influence of the choice of c σ α,0 on the V dia x drift velocity at the UOT.

FIG. 12 :FIG. 13 :
FIG.12:The influence of the grid on plasma potential at the UOT.

FIG. 14 :
FIG.14:In (a) the default grid is shown in the vicinity of the X-point.In (b) the modified grid is shown at the same location.Both grids have the same amount of grid cells, but in the modified one, no refinement is done towards the divertor targets.The rectangular B2.5 grid is shown in red, and the triangular EIRENE grid in blue.

FIG. 15 :
FIG.15:The appearing instabilities when the viscous terms are included.The instabilities on the energy crossing the separatrix for different grids are shown.

FIG. 17 :
FIG.17:The influence of averaging on the plasma potential at the UOT.

FIG. 19 :
FIG.19:The effect on the energy crossing the core boundary caused by the employed under-relaxation factors in the boundary conditions.

FIG. 20 :
FIG.20: Influence of the under-relaxation factor on j s at the UOT in attached Ne-seeded simulations

FIG. 21 :
FIG.21:The plasma potential at the UOT in detachment.

FIG. 22 :
FIG. 22: Physical mesh for the plasma in the performed SOLPS-ITER simulations.The grid boundary towards the main chamber wall is indicated in blue.Through this boundary, the energy leaving the simulation domain should be limited to draw physical conslusions from the simulation.
FIG. 23: The 2D profiles in the vicinity of the upper divertor for n e and T e in the final simulation without (a,b) and with (c,d) drifts.

TABLE I :
The number of EIRENE particles per iteration involved in the simulations.

TABLE II :
Percentage of energy crossing towards the main chamber wall towards the total energy entering the SOL.