Numerical error analysis of SOLPS-ITER simulations of EAST

Plasma edge simulations with codes like SOLPS-ITER are widely employed to interpret fusion experiments. However, numerical errors appearing in such simulations are rarely investigated, despite their potential large impact on simulation results. These errors consist of the statistical error and the bias, both resulting from the finite number of employed EIRENE Monte Carlo particles and incomplete convergence, and the discretization error due to the finite resolution of the computational grids. In this contribution, the resulting numerical errors on simulations of pure deuterium and neon seeded H-mode EAST discharges are examined. The statistical error can be kept small compared to other numerical error contributions by averaging the plasma profiles. This allows investigating the bias and discretization errors using Richardson extrapolation. It is shown that grid refinement and the number of employed Monte Carlo particles have the largest influence on the result, in agreement with similar studies of an ITER deuterium case. For the first time, numerical error bars on the entire simulated target profiles are determined showing that the largest numerical error is 17.9%, mainly due to the plasma grid discretization. On top, also numerical errors on simulated neutral pressures are investigated in detail, for which the statistical error is dominant. The analysis demonstrates which setup is needed to keep numerical errors limited: the SOLPS-ITER averaging procedure should be employed including enough EIRENE particles, and the involved grids should be sufficiently refined to reduce discretization errors.


Introduction
Current tokamak experiments focus on the use of extrinsic impurities to decrease the power and particle fluxes towards the divertor targets [1]. To investigate the effect of such impurities, a good understanding of the physics happening in the plasma edge of a tokamak is needed. Therefore, modeling is often used to examine the plasma edge in more detail [2][3][4][5].
Starting from the general Boltzmann equation, a collisional plasma can be described by the Fokker-Planck equation [6]. When the collisionality is high, a fluid approximation can be made and the behavior of the electrons and ions in a tokamak plasma can be described by the Braginskii equations, which are conservation laws for density, momentum, and ion and electron energy [7]. In these Braginskii equations, source terms appear due to interaction of plasma with neutral atoms and molecules present in the tokamak edge. As the neutrals' collisionality is not everywhere high enough to allow a fluid description, they are better described kinetically by solving the Boltzmann equation.
The B2.5 code is a finite volume (FV) code solving these Braginskii equations for the plasma particles [8]. The EIRENE code is a Monte Carlo (MC) code solving the Boltzmann equation [9] providing statistical estimates of moments of the neutral particle velocity distribution (like density). The SOLPS-ITER code package couples these two codes to model plasma and neutral behavior in the tokamak plasma edge [10,11]. While handling such numerical codes, it is important to consider possible errors on the simulation result as they can influence the solution drastically and lead to wrong physical interpretations as demonstrated in references [12,13], where for different ITER cases the choice of the numerical grid changes the simulated target profiles up to 40%. Four types of errors can be present: modeling, code, user, and numerical errors [14]. On top, also modeling uncertainties (material parameters, reaction rates. . . ) can influence the result of a simulation.
By employing SOLPS-ITER as a user, the focus is on avoiding user errors and keeping the numerical ones small. User errors appear mostly due to wrong input parameters, while numerical ones depend on the strategy adopted to solve the continuous governing model equations in a discrete setting. Moreover, the employed numerical strategy plays an important role in balancing the run-time/computational cost versus accuracy [15]. By avoiding user errors and addressing numerical ones carefully, the remaining errors present in the problem are limited to code errors, modeling errors, and uncertainties on input parameters. Uncertainty quantification through sensitivity analysis can be used to examine the impact of these uncertainties on simulation results [16]. For the examined EAST case, the relative sensitivity towards some of the unknown or uncertain input parameters has been analyzed in reference [17].
The coupled FV-MC nature of SOLPS-ITER leads to complex interactions and various error contributions, that complicate the analysis compared to standard error estimation for FV codes. Only recently these error contributions have been unraveled. Ghoos et al developed a strategy to compute the numerical errors, assessing a trade-off between accuracy versus computational time [15]. This strategy can estimate and reduce the statistical error and bias allowing investigating the discretization error. For purely deuterium ITER cases [12,18] it showed that finer grids decreased the estimated discretization error at the divertor targets with a factor of around four (from a discretization error of maximum 60% to one of maximum 15%), and displayed significantly higher heat loads on the ITER divertor than predicted with coarser grids.
This strategy is now applied to EAST simulations for which modeling was started in reference [2]. Compared to reference [12], in this study we apply the error analysis strategy to an impurity case, providing in addition error estimates for key neutral quantities such as pressures. These, together with error estimates on plasma profiles, ensure correct comparison with experimental data.
In the next section, the numerical errors appearing on simulation results are described more in detail. Afterwards, the examined SOLPS-ITER setup is discussed in section 3, focusing on the parameters affecting numerical errors. Afterwards, the errors on this particular case are analyzed in section 4. The paper concludes with general recommendations for the numerical setup derived from this study.

Numerical errors in SOLPS-ITER
In SOLPS-ITER simulations the numerical error num appears due to the involved numerics in the B2.5 FV code and in the EIRENE MC code. It is composed of four parts [15,19]: d is the discretization error due to the finite resolution of the grids on which equations are solved. c is the error arising from incomplete convergence of the residuals in the FV code. b is the deterministic (or finite sampling bias) error present due to limited number of included particles in the MC code in combination with the nonlinear coupling of the MC code with the FV one. The statistical error s appears as (statistical) noise which is linked with the finite number of MC particles. This statistical noise has a probability distribution with zero mean and standard deviation dependent on the problem setup.
In the further scope of this paper, the appearing numerical error is examined in three parts: first the statistical error s , secondly the combination of the convergence error ( c ) and the finite sampling bias ( b )-together called the bias-, and finally the discretization error d . This order is predefined. In fact, the discretization error is the error due to the finite grid spacing in case of an infinite number of MC particles per iteration [15,20]. The statistical and bias error represent the influence of the finite number of particles included in an EIRENE simulation on a given grid. The deviation of the averaged solution-the solution in which the plasma quantities are averaged over 45 000 iterations-compared to the case of an infinite number of particles per iteration causes the bias error, while the random fluctuation with zero mean in comparison with this averaged solution causes the statistical error. Therefore, this random fluctuation (statistical error) should be investigated first, before analyzing the bias error.
In the next section, the simulation approach is first discussed. Subsequently, a more in depth discussion of the different error contributions is provided.

Simulation approach
Based on the work of Ghoos et al [12], a recommended simulation approach was developed for SOLPS-ITER in reference [21] and is described in this section. The presented simulations are composed of three parts: an initial transient phase, several 10 000s iterations over which the plasma state is averaged, and a final purely EIRENE simulation in the end in which the fixed averaged plasma state is employed.
In the first, transient phase, the plasma and neutral quantities are still changing each iteration and did not yet reach their converged values (which are only effected by statistical noise). As the coupling between B2.5 and EIRENE is highly non-linear, this transient phase will typically require several 10 000s or even 100 000s iterations.
Visual inspection of key plasma quantities shows when the state of the plasma is only changing due to the MC behavior of the EIRENE code. At this point, the initial transient phase is over. Key plasma quantities monitored in this phase are: temperature at the outer midplane (OMP), and upper outer and inner target (UOT and UIT) separatrix (T OMP e,sep , T UOT e,sep and T UIT e,sep ), density at the OMP, and UOT and UIT separatrix (n OMP e,sep , n UOT e,sep and n UIT e,sep ), and maximum temperature and density at the targets.
Afterwards, the simulation is in statistically steady state [21], and it is still ran for several thousands iterations. The plasma state is averaged in post-processing over these iterations resulting in a smaller statistical error. This strategy allows using less MC particles during both the transient and averaging phase. In that way it speeds up the code to reach a converged solution [18]. Using the built-in averaging scheme only plasma quantities are stored and averaged to save memory. At the final stage, an EIRENE-only iteration with increased number of MC particles is performed on the averaged plasma background to obtain estimates of the neutral parameters with small statistical error.

Statistical error
Statistical noise as defined here has a probability distribution with mean value 0 and standard deviation σ s which acts as estimate for s [19]. It is linked to the number of EIRENE particles per iteration (P), and to the number of averaging iterations (I): σ s ∝ 1 PI [15]. This means that the iterations over which the plasma state is averaged will decrease the statistical error on the plasma quantities, while the involved number of MC particles in the final, purely EIRENE iteration determines an additional statistical error on the examined neutral quantities. For the plasma quantities, the central limit theorem leads to: σ s = σ s 500 R T 500 (2) in which σ s 500 is the standard deviation calculated on the batches of 500 iterations (the averaging iterations are divided into batches, for which a size of 500 iterations is recommended), R is the number of such batches, and T 500 is the estimated correlation time between them. R T 500 can be interpreted as that the number of batches that is large enough to be considered as independent between each other. More information about these parameters and how they are calculated can be found in references [21,22].
As recommended in reference [12], the estimated σ s is multiplied by three as there might be errors in this estimate due to the approximation of T 500 and because of the stochastic behavior of the neutral particles.
In contrast to the statistical error on plasma quantities, the statistical errors on EIRENE quantities computed during the final stage do not have any link with their value on the previous time step. This means that the statistical error on the EIRENE quantities is dominated by the final, purely EIRENE simulation step, although the (small) statistical error from the averaging iterations will have an influence on the neutral source terms. This last contribution is neglected in this work.
For one single EIRENE iteration, a root mean square value is defined in the code as approximation for the standard deviation [23]: (3) in which N is the number of included particles per iteration, X i the simulated quantity for particle i, andX the averaged quantity over all particles. In the EIRENE code, the relative form of equation (3) can always be calculated [23] and acts as the estimate for the relative statistical error (in percentage): with R g (N) the MC estimate for the studied EIRENE quantity R g , the average over N MC histories for this quantity scaled with the source strength (more details about the calculation of R g (N) can be found in reference [23]).

Bias
The bias on plasma quantities is estimated using a method similar to Richardson extrapolation [24,25] by comparing three simulations on the same spatial grid in which the number of EIRENE particles is subsequently increased by the same factor n, leading to a, na and n 2 a particles with a the original number of particles [15]. Richardson extrapolation was originally developed to study the discretization error, but in reference [15] it is demonstrated by Ghoos et al that a similar procedure can also be used to estimate the bias, as this is proportional to the inverse of the involved number of EIRENE particles bias ∝ 1 P . Richardson-like extrapolation results in the following estimate of the bias ( bias ) for quantity f a in which a MC particles are involved: with p d the order of convergence equal to [25]: and to the following estimate for bias on the simulation with na MC particles: In order to use these formulas, the bias should be the dominant error. Therefore, it can only be determined after eliminating the statistical error through post processing averaging.
As recommended in reference [12], a safety factor two is applied on the estimated bias to account for the remaining statistical error on the solution and because of possibly not reaching the theoretical bias ∝ 1 P scaling. This theoretical scaling determines that the order of convergence is one while studying the bias.
As the neutral parameters are defined by a purely MC process, there is a bias only in non-linear cases, as in linear cases the cross-correlation between the MC particles is zero [15]. As will be shown later on in this work, the effect of non-linear neutral-neutral interactions in the studied simulations is negligible, and so is their bias. Nevertheless, as for the statistical error, the bias on the coupled B2.5-EIRENE iterations will influence the neutral source terms, and in that way the final simulation result for the neutrals.

Discretization error
The last error which is investigated is the discretization error, due to the discretized equations solved on a finite grid. In SOLPS-ITER two grids are involved: the B2.5 plasma grid, and the EIRENE neutral one. This means that both of them cause a discretization error. They are analyzed using Richardson extrapolation, by performing simulations on three successively refined grids. If p d is known-which is the case in SOLPS-ITER based on the employed numerical schemes-, only two simulations are required to estimate the bias.
As recommended in reference [12] by Ghoos et al and in reference [24] by Roache, a safety factor 1.25 is applied to the error estimate to take into account inaccuracies in the estimate of p d .

Setup and result of the SOLPS-ITER simulation
In the presented work, the numerical errors introduced in the previous section are analyzed for two EAST discharges which are currently investigated with the 3.0.7 master version of SOLPS-ITER [2]: a purely deuterium H-mode discharge, and one with Ne-seeding to achieve detachment through extrinsic impurity seeding. To both discharges an external heating power of 2.25 MW is applied. Through a feedback, the separatrix density at the OMP is kept constant: n OMP e,sep ≈ 0.9 × 10 19 m −3 . The gas puff injection and pumps location is indicated in figure 1. The employed boundary conditions are based on the ones used in reference [2] and can be found in detail in appendix A.
To guarantee numerical stability and accuracy, the size of neighboring plasma grid cells should not differ too much and large cells should only be present in regions of shallow gradients of the plasma quantities. This results in three main requirements for the B2.5 grid: a symmetric grid around the X-point, smooth transitions between the different grid cells, and small grid cells at the divertor targets (due to the expected steep plasma gradients at these locations). These conditions are fulfilled in the grid shown in figure 1, using 93 grid cells in the poloidal direction, and 42 in the radial one.
The EIRENE model and reactions employed in the purely deuterium simulations can be found in reference [26]. On top, some neutral-neutral interactions from the AMMONX database are included [27]. However, their effect is expected to be limited. In the previous section, it was already indicated that the bias on the neutral quantities in the final purely EIRENE steps was neglected due to the small effect of neutral-neutral interactions, which causes non-linear effects. The probability that the involved neutral-neutral reactions take place during the studied discharge is determined by the mean free path for  elastic collisions (λ el ) [23]: in which σ el is the typical collision cross section (σ el = π d 2 4 with d the diameter of the atom/molecule) and n N the gas density. Using the maximum density in the EIRENE simulation for D 2 and D, the self collision mean free path is calculated for both substances. This results in a mean free path of ∼4.61 m for the D self collision and one of ∼0.84 m for the D 2 self collision. This means that their influence on the final simulation outcome is small as the calculated mean free path is large in comparison with the size of EAST. These findings are in agreement with the JET results of reference [26] and the TCV ones of reference [28] indicating that these neutral-neutral effects only play a role for large tokamaks like ITER, where more dense neutral regions are expected.
The gas puff for both the deuterium and the neon injection is located at the outer target strike point. More details about the setup for EIRENE can be found in appendix A.
The EIRENE code employs a triangular grid. Steep gradients of neutral density are expected around the pumping (indicated in orange in figure 1) and puffing (indicated in green) locations. Therefore, the grid is refined at these locations. The initial triangles are generated with a maximal triangle side size of 8 cm. At the locations with steep gradients, this maximal triangle side size is reduced to 2 cm. A view of the resulting EIRENE grid is given in figure 1. The grid is extended below and above the figure: the top of the grid is at the upper pressure gauge which is ∼2.5 m above the vessel wall, and the bottom of the grid is at the lower pressure gauge which is ∼3.0 m above the vessel wall.
In the used setup, EIRENE employs 10 CPU's in which the workload is evenly divided between the available processors [29]. If the number of EIRENE particles is multiplied with 10, the required wall clock time nearly doubles.
The neutral particles originate from three sources: deuterium molecules puff at the OMP, volumetric recombination in the plasma, and the surface recycling at the walls of the B2.5 grid. An overview of the included number of particle histories in the setup is given in table 1.
For all presented simulations a time step of 5.0 × 10 −5 s is employed as smaller time steps lead to a solution with the same accuracy, but slow down the convergence, and larger time steps lead to numerical instabilities. 45 000 iterations were used for the averaging phase and after the final phase the following convergence criteria were checked: • Error on the global particle balance smaller than 1.5% • Error on the global energy balance smaller than 1% • Statistical error on analyzed plasma quantities (T OMP e,sep , n OMP e,sep , T UOT e,max . . . ) obtained with the techniques of the earlier introduced SOLPS-ITER averaging smaller than 0.5%

Statistical error
The statistical error on the plasma quantities can be calculated automatically leading to the errors listed in table 2.
Two key quantities are examined in this table: the electron density (n e ) and temperature (T e ). The analyzed locations are the OMP and the UOT. These locations, which are indicated in figure 1, are the most interesting to examine in simulations as experimental data are available over here making a comparison between experiments and simulations possible. At each location the maximal simulated quantity is investigated in detail as well as the quantity at the separatrix.
On top, the statistical error on the entire T e profile at the OMP is calculated. This analysis indicates that the maximum statistical error on the T e profile at the OMP remains below ∼0.6%.
Based on table 2 it is clear that a small statistical error on the plasma quantities is obtained employing the averaging method. For a small statistical error on the neutral quantities, the number of involved MC particles in the final, purely EIRENE iteration determines the solution, while the number of particles per iteration employed during the averaging does not play a significant role.
The EIRENE quantities which are analyzed in the presented error analysis are the neutral pressures. It is chosen to analyze the neutral pressures at three locations-the ones where pressure gauges are located in EAST and until where the SOLPS-ITER grid is extended: 2.5 m above the vessel (P up ), 3.0 m below the vessel (P down ), and around the OMP at the same location where the D 2 injection is located (P mid , see figure 1). The neutral pressure depends on two EIRENE quantities, namely the atom and molecule energy density, and it is approximated as 2 3 of the sum of these quantities. As a result, also the estimate for the statistical error on the neutral pressure is the sum of the estimated errors on each quantity multiplied with 2 3 (assuming that there is no correlation between them). The locations above and below the vessel where the neutral pressures are analyzed are far away from the neutral sources (puff and main divertor). This means that only few MC particles will reach these surfaces, so that the expected statistical error is large and significantly more EIRENE particles are needed to reduce this error.
It was verified whether it was sufficient to multiply the number of EIRENE particles in the final stage with a factor 100. To study this, the number of followed particles in the MC process is further multiplied with 1000, 10 000 and 100 000. The result for the studied neutral pressures in the purely deuterium simulation is shown in figure 2(a) indicating that the simulation outcome is significantly affected by using more MC particles. From the figure it is clear that at least 1000 times more EIRENE particles are required in the final iteration. In case of Ne seeding, 100 times more particles in the final iteration are sufficient as can be seen in figure 2(b). It should be noted that the required multiplication factor for the EIRENE particles cannot be compared between the two cases due to the different plasma conditions and gas injection.
The estimated statistical errors on the pressures for the purely deuterium case are given in table 3. As the active divertor is the upper one, the particle sources will be the furthest away from the lower divertor so that fewer particles will reach the bottom measuring surface, resulting in the largest statistical error. Even when using 100 000 more EIRENE particles, the estimated statistical error on the neutral pressure is 17.4%.
For a large number of particles N, the root mean square estimate for the relative standard deviation scales as σ s (N) ∼   Table 3 and figure 3 indicate that for P up and P mid at least 10 000 times more EIRENE particles are required to ensure convergence in the purely deuterium simulation. If this number of EIRENE particles is employed, the value for s · √ N has come to a constant value. For P down , however, it is not clear if the simulation with a multiplication factor of 100 000 has already a converged s · √ N. For the further analysis in this work, 10 000 more particles were used in the final iteration.

Bias
An important condition to use Richardson-like extrapolation to estimate the bias, as introduced in section 2.3, is that the three simulations required by this method are in the linear regime. This means that the simulation result is changing linearly due to changes in the examined input parameter-in this case the number of involved EIRENE particles-so that a correct estimate of the order of convergence-p d -is possible. Therefore, simulations with 10 and 100 times more EIRENE particles are employed for the purely deuterium case, as a simulation with 10 times less particles is not in the linear regime as can be seen in figures 4 and 5. In this last figure the calculated bias is given as a function of the number of EIRENE particles. As indicated in section 2.3, the bias scales with ∝ 1 P . This means that the bias should decrease with 1 P with increased number of involved particles. In figure 5 it is clear that this is the case for the simulations with 1, 10, and 100 times more EIRENE particles, but not for the one with 10 times less particles. This is also proven by the evolution of the maximum saturation current which is given on the right axis from the figure. Also a simulation with 10 times less particles in the neon seeded case is not in the linear regime. Because the additional EIRENE particles in the neon seeded study slow down the simulation, no convergence could be reached for a simulation with 100 times more EIRENE particles. Therefore, the bias for the neon-seeded case is analyzed based on the original simulation and one with 10 times more EIRENE particles. As a consequence, assump- In black ∝ 1 P is given which shows that the simulation with 10 times less EIRENE particles is not in the linear regime, but the others are and converge faster than expected from the 1 P scaling. This is also proven by the evolution of the maximum value for j sat.,UOT which is given in red on the right axis.
tions about the order of convergence are made as discussed later on.
The bias on the plasma profiles is examined at the OMP and at the divertor targets. At the OMP, the influence of a changed number of EIRENE particles is negligible on the entire profile (both in the purely deuterium case as in the case with neon impurities) as the density at the OMP separatrix is fixed for all simulations.
As indicated in figure 4 for the profile of the ion saturation current ( j s ) at the UOT in the purely deuterium case, the calculated p d does not result in a usable estimate in each simulation point: where the influence of the bias is too small, the statistical error dominates and formula (6) cannot be used. At these locations the order of convergence is set to 0.5.
Since there are only two simulations available for the estimate of the bias in the neon seeded simulation, p d is put equal to 0.5 everywhere in this case.
The resulting bias for the inner and outer target in the purely deuterium simulation, and the one for the outer target in the neon seeded case are shown in figure 6.
This resulting bias for the plasma quantities at the inner and outer target helps to determine whether enough MC particles are involved in the 45 000 averaging steps. On the figure it can be seen that the maximal bias is obtained on the density profiles. At the UOT target, the maximum density in the purely deuterium simulation is increased by 41.5% on the bias-corrected solution and at the upper UIT even by 44.7%. Therefore, ten times more EIRENE particles are included in the simulations which are used to determine the discretization error. This results in a significant decrease of the estimated bias on the maximum density for the purely deuterium simulation down to 3.17% at the outer target, and 3.08% at the UIT. However, it should be kept in mind that this increases the required computing time by 50% when keeping the number of CPUs fixed. This illustrates the difficult balance between accuracy and required computing time. Increasing the number of EIRENE particles also reduces further the statistical errors of table 2 below 0.15% for the examined quantities. The profiles of figure 6 show that the bias is small at the locations where p d could not be estimated accurately. This, in combination with the small statistical error, allows the assumption for p d which was made before.
The bias on the neutral quantities is mainly caused by the influence of the bias on the coupled iterations which determines the neutral sources. Therefore, the effect of increasing the number of EIRENE particles in the coupled iterations is shown in figure 7 for the purely deuterium case (a) and the neon seeded case (b).
The error itself is quantified for the upper pressures in table 4. For the neon seeded case, the same order of convergence as for the unseeded simulation is assumed (p d = 1.0).
In section 2 it was indicated that first the influence of the statistical error should be investigated to ensure its influence is smaller than the influence of the bias. Comparing estimates for the upper pressure statistical error from table 3 and the bias from table 4, indicates that the statistical error is higher than the bias when 10 times more EIRENE particles are included, at least for the deuterium case. That is also why the bias for middle and lower pressures cannot be estimated with the current setup.
Similar to the plasma quantities, the bias on the neutral pressures is decreased by a factor ∼10 when using 10 times more EIRENE particles and assuming p d = 1. For the neon seeded simulations, the estimated bias stays high, but from figure 7(b) it is clear that involving 10 times more EIRENE particles changes the pressure drastically. However, as long as no simulation with 100 times more EIRENE particles is Figure 6. The outcome at the inner (left) and UOT target (middle and right) for the analyzed purely deuterium (yellow) and neon seeded (black) simulation is shown including a shaded region which represents the error due to bias (including safety factor). performed, it is not possible to verify if P up with the original number of EIRENE particles is in the linear regime, and the bias estimate can be wrong. For the further analysis, the simulation with 10 times more EIRENE particles is employed.

Discretization error
Similarly to the bias, Richardson extrapolation is employed to study the effect of spatial discretization. For the discretization error, only one additional simulation (with a different grid) is involved. The numerical scheme of the SOLPS-ITER code is a hybrid between first-and second-order numerical schemes [8,29]. This results in p d between 1 and 2. Therefore, the most conservative one (p d = 1) is taken everywhere, and only two grids are required. For this extra simulation, the already existing grid from reference [2] is used. There, the B2.5 grid is refined with a factor of 1.66 in comparison with the original grid, while the EIRENE neutral grid outside the plasma domain stays similar. Inside the plasma domain, each B2.5 grid cell is divided into two EIRENE grid cells.
In order to determine the discretization error, also the bias on this further refined grid is computed first, in the same way as for the coarser grid.
The discretization error is only analyzed in detail for the purely deuterium case as the same grids are employed for the neon seeded simulations.
The combination of the discretization error and the bias on the original simulations at the UOT target is shown in figure 8 for the standard (coarser) grid and for the refined grid. As the statistical error was less than 0.15% for all analyzed parameters, this is not included in the given error bars.
Based on these simulation results, the maximum error due to a combination of the bias and discretization can be examined. As for the bias, the maximum effect is observed on the difference in peak density. The estimated total error on the coarser grid is maximum 17.9% where the maximum estimated total error on the refined grid is 5.66%. This indicates that finer grids result indeed in smaller errors for the studied EAST cases. However, by refining the grid by a factor of 1.66, it is expected that the computing time for this 2D problem will increase with a factor around 1.66 × 1.66 ≈ 2.76. This is confirmed by the simulations. Running 10 000 iterations using 10 CPU's takes ∼5 hours for the coarser grid, but ∼12 hours for the refined grid using the same time step in both simulations. On top, the error mentioned above is the maximum error. As can be seen in figure 8 the other errors (especially on the q t and T e profiles) are much smaller. Therefore, the original B2.5 grid is used for analyzing the plasma and neutral quantities.
To examine the influence of the EIRENE grid on the plasma quantities, two additional simulations were performed decreasing the maximum triangle side size by a factor two each time keeping the same number of MC particles. No change was observed on the plasma profiles at the targets and OMP indicating that the influence of the EIRENE triangle size outside the plasma domain is limited on these profiles.
The effect of finite grid spacing on the neutral pressures is shown in figure 9 in which the simulated neutral pressures are shown for a coarser EIRENE grid (∼4 times larger grid cells outside the B2.5 grid) and a more refined B2.5 grid in comparison with the original simulation result.
It is clear from the discussion before, and from figure 9, that refinement of the EIRENE grid has a negligible impact on both plasma and neutral quantities. This indicates that there are small gradients in EIRENE quantities between the neighbouring grid cells. Earlier it was demonstrated that non-linear effects are negligible for the examined cases. If they become important, a finer EIRENE grid would be required to resolve these non-linearities.
The effect of the B2.5 grid on the neutral pressures, on the other hand, is larger. For P up the difference between using the coarse and refined B2.5 grid, is 20.4% for the purely deuterium case. Following a similar procedure as for the plasma quantities would lead to an estimated total error of 63.8%. It should be mentioned that, as no averaging is applied on the neutral quantities, a larger statistical error will influence this estimate making a correct interpretation difficult. Therefore, the total error is not estimated for P mid and P down . For the refined grid, the total error on P up stays large, but decreases to 38.3%. Despite the possible inaccuracies on this estimate, this illustrates the importance of a fine B2.5 grid to examine the appearing neutral pressures in detail. Figure 8 indicates that the magnitude of the error differs significantly between the examined plasma quantities at the UOT. In this section, insight is given on how errors on state quantities translate to errors on derived quantities. All derived quantities in SOLPS-ITER are a combination of a convective and a diffusive term as below for the velocity in the poloidal direction x [8]: (9) in which b x = B x /B is the magnetic field pitch, B the magnetic field, h x = 1/ ∇x the metric coefficient in the poloidal direction, V the velocity and D na ⊥ the anomalous particle diffusion of the species.

Discussion
Some plasma quantities are mainly dominated by the convective contribution, while others are dominated by the diffusive one. This is shown in figure 10 for the simulation results of Figure 8. Outcome at the UOT for the analyzed simulation on the original grid (left) and on the refined one (right) including the numerical error. In yellow the original simulation result is shown, the bias-corrected one is indicated with the dotted blue curve. The shaded blue region indicates the numerical error due to the bias and the finite B2.5 grid on the simulated profiles. The shown outcome is the one with 10 times more EIRENE particles.
j s and q t . The ion saturation current is nearly completely convective, whereas the conductive part has the major contribution to the heat flux. When we compare this with the bias of figure 6 and the total error of figure 8, it is clear that the numerical error is larger when the quantity is dominated by the convective term. This is in agreement with the findings of reference [30] where it is demonstrated for the discretization error that the magnitude of this depends on the convective/conductive behavior of the examined quantity.
To summarize, all the different error contributions to the main plasma target quantities and examined neutral pressures are given in table 5 using the grids of figure 1, 10 times more EIRENE particles in the coupled simulation in comparison with what was proposed in section 3, and 10 000 more particles in the last, purely EIRENE iteration. As the statistical errors on the saturation current and heat flux are not calculated automatically by the averaging tool built-in in SOLPS-ITER, they are not quantified. The discretization error is only calculated in detail for the purely deuterium simulation as the same grids are used for all simulations. On its own, the discretization error on the maximum density at the outer target is larger than the total error shown in figure 8. This is because bias and discretization shift the plasma profiles in opposite directions. Figure 6, on the other hand, indicates that the bias for the Ne seeded simulation shifts the plasma profiles in the opposite direction as for the purely deuterium simulation. As table 5 also shows a larger bias for the Ne seeded simulation, the total error can be up to ∼28% for the maximum density at the outer target. The resulting errors at the UIT are smaller in comparison with the ones at the outer target.
For the neutral pressures, the estimates for the statistical error from table 3 are an upper limit for the error in the final simulation setup as simulations with 10 times more EIRENE particles in the coupled iterations are employed. Nevertheless, it is clear from the table that this error is dominated by the discretization error making it difficult to give an accurate interpretation of the simulated neutral pressures far away from the plasma where often pressure gauges are placed in fusion experiments.

Summary
In this work the importance of an accurate choice of numerical parameters in plasma edge simulations with SOLPS-ITER is demonstrated for a purely deuterium and a neon-seeded simulation of EAST. Knowledge of this error is required to correctly interpret the simulation results and to compare them to experimental data.
In SOLPS-ITER the numerical error consists of four contributions: the statistical, bias, and discretization error, and the effect of incomplete convergence. The first two are caused by the involved number of MC particles in the EIRENE code, while the discretization error is caused by the finite grid spacing. Incomplete convergence originates from residuals in the B2.5 code.
By applying the previously developed SOLPS-ITER averaging scheme, the statistical errors on the results are reduced to negligible values below 0.15% for the examined plasma quantities in the studied discharges. The performed analysis shows that at least 10 000 times more EIRENE particles are required for the final, purely EIRENE iteration compared to the number of particles employed during the averaging. Nevertheless, the statistical error on the neutral quantities remains significantly larger than on plasma quantities, especially in remote locations reached by few particle trajectories and thus providing poor statistics.
After examining the statistical error, the focus is on investigating the bias. By applying Richardson-like extrapolation, the bias is examined and the simulation setup is modified to reduce this error contribution. This leads to a maximum bias of 3.08% for the purely deuterium simulation which is observed on the peak density at the UIT. As it is demonstrated that the non-linear effects have a negligible role for the case studied here, also the bias on neutral quantities caused by the final purely EIRENE simulation is negligible. The coupled iterations on the other hand lead to a bias of 2.14% on the upper pressure.
Finally, it is shown that the finite B2.5 grid causes the largest error on the plasma profiles resulting in a total numerical error up to 17.9% for the examined simulations. As the non-linear effects are negligible in this case, the EIRENE grid has no observable influence and is best taken coarse to decrease the computational time. The employed grids also have an influence on the calculated neutral pressures, but the error on these pressures is dominated by the statistical MC noise. However, after making a decent choice for the EIRENE particles, the discretization error shows the largest error contribution which is ∼64% for the examined purely deuterium case.
The performed analysis of the numerical errors on the plasma and neutral quantities shows that a fully converged SOLPS-ITER simulation with errors on the particle and energy balance below 1.5% and statistical errors below 0.5% does not guarantee that the remaining numerical errors are small and that the simulation result can be used for a physical interpretation. Therefore, in all SOLPS-ITER simulations the appearing numerical errors should be considered before starting with a physical interpretation of the simulation result. This has to be done taking into account the computational cost. Special attention should be put on the required number of EIRENE particles to ensure a simulation is in the linear regime for the bias. The same applies for the involved B2.5 grid. If the mean free path of the neutral-neutral collisions is large in comparison to the tokamak size-which is the case for the studied EAST discharges-computational time can be saved by excluding neutral-neutral reactions and using a coarse EIRENE grid. A more in-depth overview of recommendations to limit the numerical error in future simulations can be found in appendix B.
The SOLPS-ITER simulations were performed at the Marconi supercomputer from the National Supercomputing Consortium CINECA.

Appendix A. Detailed SOLPS-ITER setup
A.1. Boundary conditions B2.5 B2.5 solves the Braginskii equations in the edge of the plasma -in the scrape off layer and a narrow region inside the core-on a rectangular grid aligned with the magnetic field. At the core boundary, the continuity boundary condition is of Dirichlet type: a constant electron density (n e = 2.8 × 10 19 m −3 ) is imposed based on the experimental profile. A constant electron and ion heat flux (P e = P i = 1.025 MW in the purely deuterium case, and P e = P i = 0.9 MW for the Ne seeded simulation) is imposed for the energy boundary condition. The energy crossing the core boundary is in agreement with the experimental input power minus the radiated power (P in − P rad ). The Bohm-Chodura condition is applied at target boundaries for all conservation equations. Details about the specific (drift-compatible) implementation for the continuity, momentum, electron energy, and ion energy equation, can be found in references [29,31]. At the boundary of the private flux zone and the radially outermost scrape-off layer boundary, a leakage boundary condition is imposed for continuity equation prescribing the radial particle flux as a fraction of the sound speed flux which is put to 1% for the studied simulation (corresponding to a leakage factor of −1.0 × 10 −2 ). A similar condition is applied for the energy conservation (leakage factor of −1.0 × 10 −2 for ions and −1.0 × 10 −4 for electrons). For the momentum, a zero gradient boundary condition is employed. The choice of these factors is based on previous simulation setups made by Dekeyser et al [3,32]. The density at the OMP separatrix is determined by a feedback loop keeping the density at n OMP e,sep = 0.9 × 10 19 m −3 which agrees with the experimental setup as indicated earlier. Making the setup drift-compatible resulted in a lower resolution grid in comparison with the one of reference [2]. As indicated in reference [3] too narrow grid cells do not work in drift cases.

A.2. EIRENE setup
The involved neutral-neutral interactions together with their reference to the AMMONX database are summarized in table 6. In the neon seeded simulations extra reactions from the ADAS [33] and AMJUEL databases are added which are summarized in table 7. 5 In the neon seeded simulations an additional constant puff of 7.2 × 10 18 Ne particles s −1 and the same injection of D 2 take place and are included in the 5 The reactions between the higher charge states of neon are covered by the B2.5 code. simulation. For both species, the gas puff is located at the outer target strike point. Chemical and physical sputtering are not taken into account.

Appendix B. Recommendations to limit the numerical error in future SOLPS-ITER studies
A big challenge for SOLPS-ITER simulations is finding the correct balance between accuracy and computational effort. An important question is which numerical error can be tolerated for the investigated problem. In case of comparison with experimental data, their accuracy should be considered, in case only the plasma quantities are examined, less attention will be put on avoiding errors on the neutral quantities. In this section, some general statements to decrease the statistical error, bias, and discretization error are formulated while keeping in mind the required computational cost. They are based on the findings in the presented simulations, but it is aimed to make them as general as possible.
As demonstrated above, and also in references [12,18], it is beneficial to use the SOLPS-ITER averaging to decrease the statistical error, or speed up the code using less EIRENE particles to obtain a similar statistical error. In the current analysis, 45 000 iteration were immediately included in the averaging phase. In the SOLPS-ITER averaging manual [21] it is recommended to use at least 30 000 iterations. From a computational viewpoint, it is better to first check after 30 000 iterations the existing statistical errors, and decide afterwards if more iteration steps are beneficial for the averaging. By deciding if the statistical errors are small enough, also the convergence criteria introduced in section 3-especially the required statistical error below 0.5%-should be considered. A small statistical error is also required for a correct investigation of the bias.
In order to guarantee correct neutral quantities, it is of crucial importance to include enough EIRENE particles in the final, purely EIRENE step which can be verified by writing out the standard deviation on the calculated quantities. While comparing simulations and experiments, it should be ensured that the standard deviation is small in comparison with the uncertainty on the measurement at the locations where a comparison is made. For the examined example, these locations are the ones of the pressure gauges. The impact of this last iteration stage on the required computational resources is limited because this final step is not an iterating process anymore. It is better to include too many than too few EIRENE particles. For the examined simulations, it was demonstrated that at least 10 000 times more EIRENE particles are required to limit the statistical error on the upper and middle simulated pressures.
For the bias, the most important is to ensure the simulation is in the linear regime. From figure 4 it is clear that too few EIRENE particles can lead to a completely different result. Therefore, it is recommended first to run a simulation with ten times less particles. If the profiles of the original simulation and the one with less EIRENE particles are changing linearly with the difference in involved particles, the linear regime for the bias is reached.
A similar effect is observed due to the finite B2.5 grid. A too coarse grid will also lead to a simulation result outside the linear regime. In the performed EAST study, two coarser grids were tested. While a 62 × 28 grid is still in the linear regime, a 31 × 14 is not anymore. However, as the largest error contribution is caused by the employed discretization, it is not recommended to use the coarsest possible grid. For EAST simulations, the employed 93 × 42 grid or a more refined one is recommended, as the maximum discretization error on the plasma quantities stays below 20%. In case the main goal of the simulation is to compare neutral quantities with experimental data or investigate them, it is shown that a more refined B2.5 grid should be considered as the B2.5 discretization error remained the largest one on the neutral quantities in the examined case.
For the EIRENE parameters, the important neutral reactions should be investigated before starting a new simulation. In case the expected mean free path of neutral-neutral collisions is large in comparison with the expected system length scale, it is better to neglect them, reducing the computational effort. Besides the studied EAST case, this is also shown in reference [26] for JET and in reference [28] for TCV. They will slow down the simulation, and require a refined EIRENE grid. When neutral-neutral collisions are important, on the other hand, the refinement of the EIRENE grid should be investigated in the same way as the one from the B2.5 grid.