Control of particle and power exhaust in pellet fuelled ITER DT scenarios employing integrated models

The integrated model JINTRAC is employed to assess the dynamic density evolution of the ITER baseline scenario when fuelled by discrete pellets. The consequences on the core confinement properties, α-particle heating due to fusion and the effect on the ITER divertor operation, taking into account the material limitations on the target heat loads, are discussed within the integrated model. Using the model one can observe that stable but cyclical operational regimes can be achieved for a pellet-fuelled ITER ELMy H-mode scenario with Q  =  10 maintaining partially detached conditions in the divertor. It is shown that the level of divertor detachment is inversely correlated with the core plasma density due to α-particle heating, and thus depends on the density evolution cycle imposed by pellet ablations. The power crossing the separatrix to be dissipated depends on the enhancement of the transport in the pedestal region being linked with the pressure gradient evolution after pellet injection. The fuelling efficacy of the deposited pellet material is strongly dependent on the E  ×  B plasmoid drift. It is concluded that integrated models like JINTRAC, if validated and supported by realistic physics constraints, may help to establish suitable control schemes of particle and power exhaust in burning ITER DT-plasma scenarios.


Introduction
The first 2D multi-fluid edge plasma modelling studies of the current reference International Thermonuclear Experimental Reactor (ITER) design were published in 2002 [1], and assumed that ITER would be operating in a high-density L-mode regime. The model for the transport assumed constant and elevated values for diffusion coefficients as well as fixed boundary conditions at the core boundary. The plasma fuelling was achieved by a gas puff from the upper port, plus an additional fixed influx from neutral beam and pellet fuelling. In these modelling activities there was a stronger focus on the divertor operational working point, and specifically it was investigated how much radiation is necessary to radiate away Nuclear Fusion Control of particle and power exhaust in pellet fuelled ITER DT scenarios employing integrated models a large fraction of P SOL (the total power entering the scrapeoff layer (SOL) from the core, P SOL ) as is required to achieve partially detached conditions in the divertor, the latter being necessary to control heat loads towards targets and plasma facing comp onents below the material limit of <10 MW m −2 . From the modelling, scaling laws for relevant operational engineering parameters could be derived as functions of P SOL and the neutral pressure in the divertor p pfr , as well as scalings for the separatrix plasma parameters and atom influx to the core.
In a subsequent paper [2] the first consistent core-edge modelling was presented for ITER. Plasma parameters like the neutral recycling flux entering the confined region (Γ core ), separatrix density (n spx ) and temperatures (T e spx , T i spx ), which had been obtained by the previous B2-EIRENE modelling, were used as boundary conditions for the plasma core model ICPS (integrated 1.5D ASTRA core plasma code). Those conditions were expressed as functions of P SOL , the averaged divertor neutral pressure p pfr , and an extra parameter f core which takes into account the fact that core fuelling (e.g. by pellets) is more efficient. In a following paper [3] a sensitivity study was pursued to derive new scaling relations for n spx , T e spx and T i spx with the result that the temperatures were rather weak functions of all the input parameters except P SOL . The latter depends on the fusion product in the deep ITER plasma core and auxiliary heating, but also on the assumptions made about the plasma transport coefficients in the edge region. The effect of the neutral transport model on the solution was also studied in detail [4,5]. In the original B2-EIRENE simulations [1] an early version of the neutral transport model was employed (EIRENE 1996). In the updated variant of the EIRENE neutral model (EIRENE 2004) more atomic and molecular processes could be included, specifically: molecule-ion elastic collisions, molecular assisted recombination (MAR), neutral-neutral collisions (neutral viscosity employing the Bhatnagar-Gross-Krook (BGK)-model approximation) and Lyman-line radiation opacity. In the papers [4,5] it was concluded that one can use the same scalings as were derived previously for the L-mode scenarios but with a correction for the neutral pressure in the divertor p PFR allowing again for a partially detached divertor regime. Additional studies [6,7] dealt with the effect of the location of the gas puff on the efficiency of plasma fuelling and of variations in the divertor dome geometry on the divertor neutral pressure. In recent years the B2-EIRENE model has been further refined, which has led Figure 1. ITER geometry as used in EDGE2D-EIRENE simulations (in its triangularised form in EIRENE). Each EDGE2D quadrangular cell of the plasma grid is split into two triangles (unshaded area) and the void area between plasma and wall assumes no plasma. Red and green surfaces pinpoint the locations of deuterium and neon gas puffs in the model, respectively. The pink surface below the D-shaped dome mimics the divertor pump using a fixed pumping albedo. to the current ITER divertor design employing the SOLPS4.3 code package. A summary is given in the paper by Kukushkin et al [8]. In 2014 a new SOLPS variant called SOLPS-ITER was released, which also features a state-of-the-art model for edge fluid drifts and currents [9,10]. All the aforementioned ITER particle and power exhaust studies had the following assumptions in common.
• It was assumed that ITER can be operated in an L-mode like regime and that static steady state scenarios are thus applicable. • Constant radial transport coefficients were applied and specifically no explicit transport barrier was assumed; no ELM (edge-localized mode) dynamics were included. • Deep core fuelling by pellets could be treated by a steadystate assumption, thus the SOL would not be affected by pellet ablation dynamics. • The ITER divertor will be operated in a partially detached regime (a stronger level of detachment at the inner target up to full detachment) and all excessive energy in the SOL can always be removed by impurity radiation.
To achieve a fusion gain of about Q = 10 in the ITER baseline scenario, an enhanced confinement regime such as the H-mode appears to be necessary to allow for an average plasma density of 8-10 · 10 19 m −3 in the core and thus to achieve sufficient α-particle heating. The operational window is constrained by an upstream n spx which saturates with increasing external gas fuelling to 3-4 · 10 19 m −3 [2,3]. Thus, a significant density gradient must be sustained in the pedestal region. With the H-mode pedestal region becoming opaque to the neutral influx, an additional fuelling by pellet injection is required.
The numerical assessment of controlling the particle and power balance for a self-sustaining plasma burn is a difficult task, since the boundary conditions for the fields of interest (core plasma performance and divertor operational regime, including plasma wall interactions) are overlapping. Top: radiated power P rad of neon-seeded impurities integrated over the entire EDGE2D-EIRENE simulation domain (SOL & edge) as a function of core-sided ion flux Γ ion (a) and neon-seeded gas flux Γ Ne (b). Corresponding peak power density q peak at the HFS/inner target (middle, (c) and (d)) and LFS/outer target (bottom, (e) and ( f )). Transport assumptions: blue-moderate confinement; red-good confinement.
Specifically, the modelling of external fuelling mechanisms by gas puffing or pellet ablation does not allow a clear separation of the plasma core, edge, and SOL. Hence, an integrated modelling approach is recommended, allowing for a consistent communication of the plasma conditions on closed and open magnetic surfaces at the same time.
In a previous study of the pellet fuelling mechanism, the advantage of an integrated simulation scheme was demonstrated [11]. When using core plasma codes like JETTO or ASTRA one normally assumes fixed boundary conditions at the separatrix, albeit pellet ablation is a transient process with finite timescales and the separatrix conditions must evolve in time accordingly.
2D SOL plasma simulations also usually assume fixed boundary conditions either at the separatrix or a little deeper in the core (up to the top of the pedestal). However, artificially fixing power and particle fluxes at the core sided boundary in edge plasma simulations ignores the transients at the time of discrete pellet ablations (and ELMs in case of H-mode scenarios) and thus prevents any assessment of the dynamical response of the SOL and divertor plasma on the confined plasma and vice versa. In this paper it is shown that this dynamic evolution of upstream parameters in the modelling of power and particle exhaust for ITER cannot be ignored.
In a complementary work to this paper, the ITER baseline scenario was modelled using a core plasma model only, which neglected the evolution of the scrape-off plasma layer at the time of pellet ablations and also excluded the divertor recycling process [12]. The model also assumed fixed boundary conditions at the separatrix for neutral influx and n e sep . In [12] it was conjectured that in case of a relaxation of the separatrix condition constraints, the dynamics of separatrix density and temperature should strongly affect the recycling process and the particle content in the SOL when pellets are ablated in the plasma core.
We will present in this paper the extension of the JINTRAC integrated model [13,14] that can be more helpful in understanding the evolution of a pellet fuelled ITER core plasma in the nuclear phase by relaxing the boundary conditions, employing a realistic 2D edge plasma model that is consistent with the previous ITER SOLPS divertor design studies. The full JINTRAC model has been successfully used to model the dynamics of JET ELMy H-mode scenarios, including pedestal fuelling (see [15,16] and references therein). Here the results of an assessment of the impact on the ITER divertor conditions on variations of P SOL due to pellet ablation and α-particle heating are presented. It will be shown whether those conditions can be made compatible with divertor material limits, i.e. whether a regime can be found in which the core confinement is good enough to achieve Q = 10 and, at the same time, the target heat loads do not exceed the critical limit of 10 MW m −2 employing a partially detached divertor regime.

Summary of sections
In section 2, the setup of the 2D SOL/divertor model EDGE2D-EIRENE, which is a central part of JINTRAC, is described for the investigated baseline H-mode ITER scenario. The underlying assumptions for the simulations and model parameters are explained. (The results of a benchmark between the used EDGE2D-EIRENE and B2-EIRENE/ SOLPS4 models simulation for an ITER L-mode case are presented in the appendix.) In section 3 the results of initial parameter scans performed with EDGE2D-EIRENE are presented, highlighting the parametric range of particle throughputs for main-plasma and neon species that are required for significant SOL power loss to achieve partial detachment. These results are taken as an additional constraint into the JINTRAC model. Section 4, the main part of this paper, presents the results of the fully integrated model, highlighting the dynamic evolution of the plasma properties caused by perturbations due to the pellet injection and α-heating and their impact on the SOL and divertor conditions. Section 5 concludes the paper.

Setup and benchmark of the 2D SOL/divertor numerical model for ITER
The SOL/divertor model EDGE2D-EIRENE [17][18][19] solves the 2D multi-fluid Braginskii equations for particle, parallel ion moment (for all ionic species), as well as internal energy equations for electrons and ions. It is assumed that the magn etic equilibrium is fixed at all times, i.e. the same 2D topology (equilibrium, vessel geometry) as used in the previous work done with SOLPS4 [4,5]. Figure 1 displays the computational mesh. The cells are aligned along the magnetic field lines in the 2D poloidal plane. By doing so, parallel and perpendicular directions are separated and thus the transport models for each direction can be defined independently from each other. It has been estimated on the basis of simple ion-gyro radius scaling that the pedestal depth (i.e. the edge transport barrier (ETB) width) at the outer mid-plane has an extent of roughly 6 cm, and thus the magnetic grid just extends to that limit in the EDGE2D-EIRENE model. All conditions for the inner core boundary apply, such that the innermost ring refers to the top-of-pedestal location. The strike points are in contact with the vertical target structures, as seen from figure 1.
To provide source-and sink-terms for the particle, momentum and temperature equations in EDGE2D, EIRENE is used to model the neutral particle distributions (for the main species and impurities) self-consistently by the kinetic Monte Carlo method [20]. EIRENE as it is used here follows test particles by a bookkeeping of successive events (e.g. recycling or collisions of neutral particles, as well as ion-neutral interaction like ionisation or charge-exchange) on a triangular 2D mesh of the poloidal cross-section. The EDGE2D plasma grid is triangularised and completed with additional triangular cells between plasma grid and vessel structure to allow for neutral-wall interactions (see figure 1). This method preserves geometric details of the ITER design which are important in estimating the neutral pressure in the divertor [8]. For technical reasons the current model does not allow the plasma to be in contact with the outer wall, and it is assumed that the radial plasma flux crossing the outermost boundary of the plasma grid is instantaneously transported to the wall and recycled there.
The parallel transport in EDGE2D is assumed to be of thr classical Spitzer-Härm kind [21]. Thus, heat conductivities κ e, i and viscosity η i along the field are strongly dependent on the local temperature (~T 5/2 ). Strictly speaking, the validity of the classical approximation for parallel transport can be expected only in cases where SOL plasma collisionality is high enough. For situations with very high temperatures and/or strong parallel gradients, kinetic corrections can be applied. In this case, so-called heat-flux limiting factors can be imposed [22], for example on parallel heat flux q || . In order to be in line with previous SOLPS4 design studies, EDGE2D was set to have a flux limiting factor for q e|| of 0.3, whilst no limit for q i|| was imposed. For the parallel ion momentum balance a viscosity limiting factor of 0.5 was applied.
Anomalous turbulent-advective transport in a perpendicular (or radial) direction is prescribed by a set of radially varying particle diffusivities D perp , viscosities η perp and (pinch) velocities V pinch defined for each ion species, as well as temperature diffusivities χ e and χ i . The actual values for these quantities (and spatial dependence) are usually regarded as free parameters. Depending on the scenario confinement assumptions (e.g. L-mode or H-mode) the absolute values of those coefficients are adjusted accordingly. The impurity particle transport model (e.g. neon) is assumed to be the same as for the main  deuterium plasma particles. The thermal force on impurity ions is included in the ion momentum balance. Energy losses by line radiation are included. As with the previous SOLPS4 studies, the effect of fluid drifts is not included in this paper.
The combined set of multi-fluid equations for the main plasma and impurities is accompanied by a set of prescribed boundary conditions. In the perpendicular outward direction these have to be defined at the outer SOL and private flux zone boundary facing the main chamber and divertor, and radial decay lengths for density and temperature are given (3 cm). The resulting plasma flows crossing the radial boundary are recycled at the vessel wall, as described above. In the parallel direction within the SOL the plasma is in direct contact with the target plates. It is assumed that the plasma flow is 100% recycled as neutral particles. Constant sheath heat transmission factors for ions and electrons are assumed here (γ i = 2.5, γ e = 4.5 using the definition as given in [22]). Secondary electron emissions are neglected. The Bohm criterion for the parallel velocity (M ⩾ 1) is assumed.
Unless the SOL solution is not coupled to a core model as described in section 4, the power and particle fluxes of electrons and ions at the inner core boundary (i.e. the innermost magnetic surface covered by the fluid mesh in the simulation) entering the edge region of the plasma P edge must be prescribed (imposing auxiliary and fusion power as well as radiation loss). In such conditions, the particle flux at the inner core boundary is also a combination of particle flux from the neutral beams and (constant) pellet-fuelling contributions, plus any additional particle ionic return flux of atoms which reached the inner boundary.  A density scan can be pursued by varying the molecular deuterium gas puffing rate applied at a surface located at the top plane (see figure 1) in order to steer the particle throughput. Additionally, seeded impurities (neon) are introduced as atoms (e.g. from the divertor region), again with a prescribed seeding rate. A pumping surface is located below the D-shaped dome at the bottom plate in the divertor (see figure 1). This toroidally continued surface is specified such that with its given extent in the poloidal plane the neutral particles have a certain probability of being absorbed. A fixed albedo coefficient, ζ albedo = 0.94, is imposed. The actual pumping speed is a function of neutral pressure and should be regarded as a dynamic quantity when large ELMs or discrete pellet ablations occur (see section 4). In order to be consistent with the previous modelling activity done with SOLPS4 [1][2][3][4][5][6][7] a dedicated benchmark of the EDGE2D-EIRENE transport model was pursued. The results are presented in the appendix. A more precise benchmark of EDGE2D-EIRENE and SOLPS4 for ITER SOL/divertor plasmas has been elaborated recently [23].

Constraining the JINTRAC SOL/divertor model: power dissipation for ITER scenarios using neon seeding
The design of the full JINTRAC model setup should focus on the dynamics of the plasma core performance (density evolution and α-heating) and at the same time on the impact and response of the (detached) divertor conditions. This numerical task is large and complex, and reasonable simplifications may be helpful. For instance, the full picture of impurity transport is not necessarily required as long as one can impose some plausible constraints on the power loss function at the edge.
Assuming that a simple temperature dependent corona model for the cooling function is sufficient (see next section), some specific questions arise. What is the total radiation loss P rad required in ITER to ensure (partial) detachment in the divertor, i.e. to keep target power densities below the material limit of 10 MW m −2 ? For a given set of input power P edge , what are the ranges of particle throughputs of neon impurity particles to be imposed in the model? (This question implicitly assumes that neon has no impact at all on the plasma core.) To provide answers, the EDGE2D-EIRENE model as described in the previous section was employed for the steadystate conditions of an ITER case. The transport coefficients were adjusted to be consistent with the assumption of the ITER core plasma baseline performance. From the previous core-only modelling [26] it was possible to derive two suitable sets of ETB transport coefficients (assuming a 'quiet' H-mode scenario in steady state with the so-called 'continuous ELMs' model, as explained in the next section): These sets of transport parameters were applied only in the ETB of a given width (6 cm at the outer mid-plane) with a narrow prolongation into the SOL near the separatrix (another 0.5 cm at the outer mid-plane), to match the typical H-mode heat flux tube width upstream as proposed in [24,25]. In the far-SOL the radial transport was assumed to be significantly higher, i.e. of turbulent or Bohm type: . not taking into account any dependence on local plasma parameters such as SOL collisionality [27] or any other spatial variation towards the far-SOL). The radial viscosity η SOL was assumed to be proportional to particle diffusivity. Finally, all transport coefficients defined at the outer mid-plane were mapped along the poloidal direction, along the flux surfaces in 2D. It was assumed that the power crossing the inner core boundary was constant in time and in total P edge = 80 MW (shared between electrons and ions at a ratio of 1: 2 respectively) as proposed in [26] for ITER scenario-2. For the variation of the particle fluxes, two scans were performed: • scan 1: a core ion-flux scan between Γ ion = 2 · 10 21 s −1 …1 · 10 23 s −1 where the puffed neon gas flux was fixed: Γ Ne = 10 19 s −1 ; • scan 2: a neon gas flux scan between Γ Ne = 1 · 10 19 s −1 …8 · 10 19 s −1 where the core boundary ion-flux was fixed: Γ ion = 2 · 10 22 s −1 .
The first scan helps to estimate the sensitivity of a varying core particle flux on detachment (i.e. from a constant flux from the pellet-ablated material in the core). The second scan is designed to find a suitable regime of radiation losses P rad by neon in the SOL at a given D particle throughput, to reduce significantly the power density at the targets. An additional molecular deuterium gas flux was assumed always to maintain a minimum level of neutral pressure in the divertor (Γ gas = 1.4 · 10 23 s −1 , i.e. the value taken from the previous EDGE2D-EIRENE-SOLPS4 benchmark exercise; see the appendix). Figure 2 shows the variation of the separatrix density and temperature as a function of the effective ion flux Γ ion from the core (scan 1). Separatrix density n e spx does increase moderately with Γ ion whereas T e spx is barely impacted by the core flux. However, whilst there is not much impact on n e spx from the change of transport (from moderate to good confinement), T e spx is strongly affected and drops by 30-40%, irrespective of the value of Γ ion .
Within the divertor and with an increasing seeding rate of neon impurities Γ Ne , the radiative power losses P rad increase significantly, as seen from figure 3 (P rad integrated over the whole simulation volume). On the other hand, increasing the value of Γ ion while keeping Γ Ne fixed does not change P rad at all. Transport also plays no significant role in P rad . Figure 3 also displays the resulting peak heat flux density q peak on the inner and outer target plates (for both scan 1 and scan 2). A first observation is that q peak is strongly dependent on the assumed radial transport model. At low Γ Ne and in the good confinement case, the outer target receives maximum power densities far above the critical limit of 10 MW m −2 and the inner target stays below the limit only if Γ ion is larger than 6 · 10 22 s −1 . The situation is better in the moderate confinement case, where both inner and outer targets stay just below the critical material limit for all values of Γ ion . Introducing neon impurities effectively dissipates excess energy from the SOL by radiation. For a given scenario of H-mode confinement, the neon gas flux can almost always be chosen such that q peak drops below the critical material limit of 10 MW m −2 . In combination with elevated values of Γ ion , even in the good confinement regime q peak can be managed.
All results presented in this section have assumed constant Γ ion and P edge . For the more complex situation that includes dynamic response to pellet ablations and transients in αparticle heating and transport, we may now assume that a minimum P rad > 40 MW is required to keep the outer divertor at least partially detached.

Integrated modelling of ITER scenarios with JINTRAC
Within the JET integrated transport code suite JINTRAC [13,14], the 1D core code JETTO [28] and the 2D SOL code EDGE2D-EIRENE [17][18][19][20] can be run in a combined way (the so-called COCONUT operation mode [29]) by interfacing the two codes via a common boundary surface which is assumed to be the separatrix. Heat and particle fluxes and/or plasma parameters can be exchanged at this communicating surface by calling JETTO and EDGE2D-EIRENE in an alternating manner. A self-consistent time-dependent solution for both core and SOL plasma can be achieved by this method. The use of such a method does not necessarily mean that a solution for the posed numerical problem can be found more easily, when compared to the problem of finding simple scaling laws as described in [2][3][4][5]. Rather, the problem is relocated: the finding of a good set of boundary conditions for each separated code has now been replaced by the finding of a suitable set of plasma models (including their transport characteristics and their own set of boundary conditions) for both SOL and core at the same time. A clear advantage of such a methodology is that the inclusion of dynamic events like ELMs and pellet ablation in a timewise evolution of the plasma can be implemented more straightforwardly.
Before presenting the results of the integrated model for an ITER H-mode scenario with pellet ablation dynamics  Time transients of critical divertor detachment characteristics for the modelled ITER high-density baseline scenarios case with P rad = 60 MW radiative fraction in the SOL and top pedestal n ped = 1.05 · 10 20 m −3 . From top to bottom: (a) relative pressure ratio f p = p omp / p targ mapped at outer mid-plane (separatrix located at R ≈ 8.19 m), (b) log T e (eV) at targets (along S coordinate in (m)), (c) target ion-recycling flux Γ targ (10 24 s −1 ) and (d) total heat flux density arriving at target plates q targ (MW m −2 ) (for (b)-(d) the separatrix is located at S = 0). The time-axis starts at the beginning of the JINTRAC simulation (t = 0 s) which began from a non-detached plasmastate. As simulation time advances, the divertor plasma is driven into a detached state (as seen from the particle flux Γ targ roll-over and drop in q targ ). In the saturated phase after 0.5 s, regular pellet injections are triggered when the main plasma falls below a minimum density n ped = 1.05 · 10 20 m −3 . These cause a transient reduction of the degree of detachment shortly after pellet ablation times.
included, we will describe briefly the setup of the coupled JINTRAC transport model.

JINTRAC transport model setup
The core fluid equations are solved by JETTO on a 1D discretised grid up to the separatrix. The JETTO's internal equilibrium solver (ESCO) [28] is used to model self-consistently the magnetic field with the evolving profile for the current.
For the anomalous transport in JETTO the empirical Bohm/gyroBohm transport model adapted to ITER conditions as described in [12] is used (using similar assumptions for the model's shape factors: linear in the core and sublinear at the edge). The neoclassical part is treated by using the NCLASS package [30]. An ETB is imposed, having a reference extent of 6 cm at the outer mid-plane. The actual level of transport at the edge is given by applying a small multiplicative factor of 0.001 to the Bohm/gyroBohm transport coefficients D anom , χ e, anom and χ i, anom to suppress the anomalous part in the ETB, in effect leaving only the neo-classical part within the 6 cm extent.
As shown in the previous work of modelling the ITER core H-mode plasma with JETTO, CRONOS and ASTRA [26], the ELMy H-mode scenario can be implemented by enhancing the transport at the edge, using an effective increase of transport coefficients to include magnetohydrodynamic (MHD) destabilising effects by ELMs in a time-averaged way (the "continuous ELM model" as described in [12,26]). To achieve this, a maximum allowable pres sure gradient parameter α c within the ETB is imposed to average over many ELM events and the transport coefficients are self-adjusted or clamped by the following expression (e.g. for the particle diffusivity): where α ~ q 2 /B 2 · ∇p is the normalised pressure gradient at the edge [12], c is a constant driving the transport if α > α c , and b is a shape parameter expressed by an exponent (here, β = 2.5). A value of α c = 1.7 is taken, this being consistent with the previous core-only predictions from [26]. Other MHD effects like sawteeth are precluded from the present analysis. In the deep core the particle diffusivities follow the harmonic average of the heat conductivities according to the Bohm/gyroBohm prediction [31].
For the ITER reference scenario a constant auxiliary power of 33 MW is applied by using the neutral beam model PENCIL [32] which is part of JETTO. No radiofrequency (RF) heating is assumed in the model. The fusion power is given by the SIMOD fusion product model [33].
The core plasma transport model deals only with fuel species (i.e. D and T ). Impurities are modelled only interpretatively by assuming a radially constant Z eff = 1.7. A radially constant power loss profile due to radiation is thus imposed, normalised to a fixed value of ≈ 43 MW in the core.
In the coupling to EDGE2D-EIRENE, the radial transport near the separatrix is assumed to be dependent on the level of transport provided by JETTO at the separatrix location. In many edge modelling attempts (e.g. at JET [24,25]) a prolongation of the ETB into the SOL was found to be required when the experimental radial profiles (i.e. the power flux width in the SOL) of the inter-ELM phase of an ELMy H-mode discharge were to be recovered by an edge model. Similarly, for ITER we assume an extra transport barrier width in the SOL near the separatrix of 0.5 cm extent and near the outer mid-plane where anomalous transport is suppressed too and is prescribed by the lower (neo-classical) level of transport provided by JETTO. Further out into the far SOL the transport is increased to some Bohm level of anomalous transport: D far-SOL = 0.3 m 2 s −1 , χ i far-SOL = χ e far-SOL = 1 m 2 s −1 . Parallel transport in EDGE2D is again assumed to be classical, as before. To incorporate kinetic effects, electron heat flux-limits are again used in the simulations. The assumed sheath heat transmission coefficients at the targets are 2.5 and 4.5 for ions and electrons respectively, as described in the previous section.
EDGE2D at its present state is not able to treat more than one fuel species at a time (independently from additional intrinsic or extrinsic impurities) and thus cannot model D and T simultaneously. However, JETTO is providing EDGE2D with D and T ion particle outfluxes at the separatrix. In the coupling procedure EDGE2D, both D and T fluxes are combined at the separatrix into a single D-flux (Γ D In turn, EDGE2D-EIRENE passes back the neutral deuterium recycling flux crossing the separatrix and D ion densities to JETTO, and a 50/50 split into D and T species is assumed in the back-coupling.
The particle, momentum and energy sources and sinks due to neutral and plasma-wall interaction in EDGE2D are derived by the kinetic Monte Carlo code EIRENE, as before. A cold neutral molecular gas-puff is applied at the top plane of the ITER vessel. The gas flux Γ D2 is assumed to match the constraints derived earlier in the benchmark with SOLPS4, as described in section 3 and in the appendix. The neutral flux eventually crossing the separatrix into the core is taken over by the 1D FRANTIC code, which approximates the particle sources from neutrals in JETTO by using a reduced diffusive fluid model.

Time-dependent modelling of discrete pellet ablation
In order to perform integrated simulations taking into account the transient evolution of plasma parameters during pellet ablation and subsequent density relaxation processes, JINTRAC has been coupled with the hydrogen pellet injection code HPI2 [34,35], which provides the pellet source profiles for a given pellet injection configuration and plasma state. HPI2 determines the change in plasma density and temperature after pellet injection, taking into account a pellet ablation model which is based on the neutral gas and plasmoid shielding (NGPS) description (allowing for electrostatic sheath effects, the spectral energy variation of the incident particles, and the time evolution of the plasmoid from its formation until its release), and a four-fluids Lagrangian model for the homogenisation process of the released drifting plasmoid clouds, including drift damping effects related to Alfvén wave dissipation and both internal and external short-circuiting of parallel currents. Following the pellet along its injection path, the ablation and homogenisation processes are calculated for a representative subset of plasmoids during their lifetimes to determine the local pellet ablation rate and drift displacement.
In an iterative manner, the locally ablated and drifted pellet particles are deposited in the background plasma, and the  plasma density and temperature profiles are updated accordingly. The designated ITER pellet launcher configuration for injections from the high field side of the plasma [36] was used to inject pellets with a mass of ≈ 6 · 10 21 atoms (50% D, 50% T ) at a speed of 300 m s −1 . The effect of plasma precooling due to previously deposited pellet particles in the ablation process leading to a reduction of the ablation rate was taken into account, although it is considered to be weak, as the pellet penetration is expected to be small with respect to pellet particle penetration after homogenisation of the plasmoids. Either 50% (as a conservative estimate) or 100% of the predicted plasmoid drift displacement were considered in the calcul ation of the pellet source profiles. As noted in section 3, impurity transport was precluded in the full JINTRAC simulation. However, additional radiation in the SOL caused by impurities was taken into account in EDGE2D by a prescription of either P rad = 40 MW or 60 MW of additional radiation losses, in order to maintain partially detached divertor conditions. The spatial radiation loss distribution was governed by a simple corona model (Neuhauser-Garching model) for each grid cell k. For a given electron density and temperature in a cell k we define: k rad e 2 e 3 2 e 3 / and for each cell k, each S k rad is rescaled to match the specified integral impurity radiation loss P rad : where V k is the volume of a cell k.

Results
Three density levels of ≈ 8, 10, and 12 · 10 -19 m −3 on top of the pedestal were targeted in the simulations (see figures 4-6 for cases with P rad = 60 MW of radiated power in the SOL, and figures 7 and 8 for cases with P rad = 40 MW). The pellet injection frequency is determined by a density feedback control mechanism which tries to maintain the density on top of the pedestal above a specified level, to ensure that the plasma particle content averaged over a pellet injection cycle remains roughly constant. The resulting pellet injection frequency ranges between ≈ 3 and 6 Hz (which is slightly higher than the results in [12], where a 25% larger pellet mass was assumed). Clearly, the pellet frequency must be raised to reach higher core densities, as particle transport is enhanced with improved fuelling. This is a consequence of the slightly increased density gradients for high density regimes, and, more importantly, a rise in fusion power leading to higher heat fluxes, which for similar temperature gradients are associated with increased heat conductivities, and, with the Bohm/gyroBohm model assumption /( ), also with larger anomalous particle transport.
For the medium density case with 50% of plasmoid E × B drift velocity and P rad = 60 MW of radiation losses in the SOL, a quasi-stationary regime could be obtained with H 98y ≈ 0.95, P fus ≈ 450 MW, Q ≈ 13, and W th ≈ 330 MJ (see figure 8), whereas for the 50% drift/40 MW SOL radiation cases, quasisteady state conditions were found at lower density with H 98y ≈ 1.0, P fus ≈ 420 MW, Q ≈ 12.5, and W th ≈ 310 MJ (see figure 11). With stiff transport models such as GLF-23 [37], the predictions for the core confinement with consideration of pellet fuelling are slightly less optimistic [26], and the energy content and fusion power are found to remain almost constant within the range of particle content that has been scanned in this study, due to pressure profile stiffness.
Regarding the need to mitigate the heat flux arriving at the target plates, the simulation cases that are analysed here are very demanding, as the predicted plasma performance translates into high total power fluxes through the separatrix, which varies transiently in the ranges ≈ 60-120 MW and ≈ 55-110 MW with 60 MW and 40 MW, respectively, of SOL radiation and 50% of the predicted plasmoid drift for all density levels. Electrons and ions carry approximately the same fraction of the total heat through the separatrix (the ratio between electron and ion heat fluxes remains within an interval of ≈ 0.8-1.1). The time-averaged total target heat flux is in the ranges of ≈ 70-80 MW and ≈ 50-70 MW with 60 MW and 40 MW respectively of SOL radiation. For the two aforementioned quasi-steady state cases, the plasma core is mainly heated by ≈ 85-90 MW of alpha thermalisation power, 33 MW of NBI power, and only ≈ 2 MW of ohmic power. It is cooled by ≈ 43 MW of core radiated power, but also, due to the interaction with incoming neutral particles in the vicinity of the separatrix, causing an extra net power loss of ≈ 10-20 MW.
The transport after pellet injection is increased by the pellet-induced perturbation of the plasma background profiles, leading to a variation in the total power flux crossing the separatrix during a pellet cycle of up to ≈ 40-50 MW. Both χ e, i and D are found to increase temporarily by up to 50% after the injection of a pellet in the outer core region that is affected by pellet particle deposition and also in the ETB zone, the latter reflecting a pellet-triggered transient increase in the ELM frequency; or, in more extreme cases, the appearance of a burst of ELMs, as was observed in the experiment [38]. The predicted pellet-induced heat flux variation may be underestimated, as it was found in [31] that particle diffusivity increases significantly above the level predicted by the Bohm/ gyroBohm model in the post-pellet phase in JET experiments. Such an effect has not been considered in these simulations. The fusion power fluctuates by only a few percent due to pellet injection, as the fusion process is dominant in the plasma core, whereas the pellet perturbation affects only the plasma edge region. The variation in P fus does not influence the heat flux to the SOL, as the change in α-heating is damped by the time required for the α-particles to thermalise. In addition, the heat fluctuation is further attenuated when it diffuses towards the separatrix.
Comparing the simulation cases at 60 MW with those at 40 MW of radiated power in the SOL at the same density (and 50% of the predicted plasmoid drift), one finds that the plasma confinement is slightly improved for the latter. This seems to be a direct consequence of the variation in SOL conditions. At lower radiated power, the upstream temperatures of both the electrons and the ions are increased. The electron temperature is directly affected by the reduced radiation level, whereas the higher energy of the ions seems to be mainly related to a reduction in heat loss via charge exchange interactions with cold neutrals, as the target plates are less well detached due to the higher incident heat flux, causing a reduction in neutral density. Whereas the electron density at the separatrix varies within a range of ≈ 4.5 ± 1.5 · 10 19 m −3 , the electron and ion temperatures are increased by ≈ 35%, from ≈ 1.0 ± 0.4 keV and ≈ 0.9 ± 0.3 keV, respectively, to ≈ 1.35 ± 0.25 keV and ≈ 1.22 ± 0.12 keV, respectively. (It is recalled that the ETB is assumed to be extended into the SOL, i.e. these values do not correspond to the temperatures at the outer edge of the ETB.) With the model constraint that the critical normalised gradient in the ETB remains unchanged, and due to the low stiffness in the predictions of the Bohm/gyroBohm model combined with the temperature sensitivity of P α , this increase in temperature can spread towards the plasma core and yield a slightly higher energy content.

Sensitivity on pellet ablation plasmoid drift
As the plasmoid E × B drift displacement is expected to be in the order of several tens of centimetres as a consequence of high initial plasmoid pressures related to strong ablation rates in a high energy plasma such as ITER, the pellet particles are deposited comparably deep in the core plasma with a deposition barycentre of ρ ≈ 0.8 for 100% of the predicted plasmoid drift, leading to improved particle accumulation in the core. As the fuelling efficiency is improved for deeper penetration of the pellet particles, the plasma density is increased in simulations where 100% of the predicted E × B drift is taken into account as compared to those considering only 50% of the drift displacement. This leads to an amplification of particle transport caused by higher density gradients and fusion power, as explained above, and an increased pellet particle throughput requiring a higher pellet injection frequency to maintain the plasma particle content. The value of H 98y remains close to unity, because the heating power and the thermal energy content increase at the same rate.
All results presented so far (figures 4-8) show time traces and profiles for simulations with the three density levels for the 50% E × B drift case. Figures 9 and 10 compare the medium density cases with 50% and 100% of E × B drift at a radiated power of 60 MW and 40 MW in the SOL. As expected in the cases with 100% of E × B plasmoid drift, the density is higher and thus Q increases by at least 30%, whereas the confinement in terms of H 89y stays more or less the same.

Divertor operation compatibility and plasma detachment stability
In the model, a discrete pellet ablation and particle deposition event has two important consequences for the SOL: (a) an increased particle flux into the SOL α SOL due to the transient increase in density and transport in the plasma edge region that is driven by the deposition of pellet particles, and (b) an increase in the ion and electron power entering the SOL (P SOL , see figures 5 and 8) due to the pellet-induced perturbation of profiles in the same region. It is stressed again that the increase in P SOL is not due to the increase of the α-particle heating P fus during pellet fuelling itself (as the α-particles thermalize over longer timescales), but due to a transient increase of convective heat flux with a localised increase of the edge density. The duration of this transient increase in P SOL depends also on the selection of parameters c and b in the continuous ELM model (equation (1)) clamping the profiles to a prescribed pressure gradient α c . The impact of transients in P SOL and Γ SOL on the divertor conditions is described in the following.
We define the following as decisive parameters in characterising plasma detachment: (a) the total pressure (static + dynamic) in the SOL shall drop significantly along the field lines, i.e. the total pressure ratio upstream/downstream should drop below unity: f p = p target / p omp < 1; (b) in front of the target there shall be a volume recombination zone which necessitates an electron target temperature of less than 2-3 eV, and; (c) the particle flux arriving at the target plates shall drop significantly towards zero (a flux roll-over). Figure 11 shows the time transients of all three of those parameters at the inner and outer target plates for the highdensity case with 50% plasmoid drift and fixed radiated power in the SOL of P rad = 60 MW. It is important to note that at the beginning of the plotted time transients (and the consecutive plots discussed below) the density is still rising up to satur ation, where the density control mechanism applied results in a regular pellet frequency. In this manner the particle flux roll-over can be clearly seen in the starting phase. The electron temperature drop is also clearly seen. Both targets begin to reattach at the time when a pellet is ablated. In the integrated model a strong increase in the power entering the SOL P SOL coincides with a strong increase in edge heat transport (see figure 5) and thus leads to a strong decrease of the level of detachment at pellet ablation time, which is ultimately responsible for the approach to the critical material limit of 10 MW m −2 at the outer target plate (see table 1).
For the medium density case with 60 MW SOL radiation this effect is somewhat moderated, as shown in figures 5 and 12. The inner target stays detached whilst the outer target reattaches only at pellet ablation time. This behaviour (although being good from the viewpoint of divertor particle and heat control) is apparently unexpected since one would assume that with decreasing density the overall level of detachment would degrade accordingly. This is contrary to conventional no-α-heating tokamak simulations where the effect is different, assuming a fixed heating power and constant P SOL ; with increasing density the level of detachment would increase further, up to a density limit. However, with α-particle heating from the fusion process in place, P SOL is positively correlated with the evolution of density in time (and on top affected by transients due to pellet ablation) and thus we observe the effect that the level of detachment actually decreases with core density.
If one goes to very low core/pedestal densities, the fusion product drops significantly. Here, no more pellets are added to the system due to the fact that the pedestal density is rising monotonically (in this case the density control mechanism will trigger a pellet only when n ped ⩽ 0.7 · 10 20 m −3 ) as seen from figure 4, and P SOL becomes saturated. The latter is such that it will not heat up the SOL strongly enough, and due to the elevated particle content in the SOL the pressure balance along the field lines leads to both targets being completely detached, as seen in figure 13. One could regard this lowest core/pedestal density scenario as uncontrollable since both targets are completely detached and already beyond its system-dependent density limit.
If the radiation in the SOL is assumed to be lower (P rad = 40 MW instead of 60 MW), the heat flux that reaches the target plates is significantly increased, for obvious reasons. In all the cases with a low SOL radiation of 40 MW that have been modelled, the maximum local heat flux arriving at the outer target in the vicinity of the line of intersection with the separatrix transiently exceeds by a factor of 2-3 the critical limit of 10 MW m −2 during or immediately after pellet ablation (see table 1). A trade-off must therefore be sought for the optimum adjustment of P rad , between the need to protect the target in the post-pellet phase and the necessity to avoid complete detachment between pellet injections. The optimum P rad depends on numerous aspects that have not been properly considered in this study, such as impurity transport, the trigger conditions and discrete nature of MHD instabilities like ELMs, and the impact of the latter on the level and distribution of divertor neutral pressure and pumping. It may thus deviate considerably from the above-cited estimate, although the inferred qualitative behaviour and dependencies of the divertor operation conditions can be expected to be valid.

Cyclical operation regimes and particle throughput control
In figure 14 n e, spx and T e, spx are displayed as a function of the instantaneous total ion flux Γ perp from the core crossing the separatrix for both the P rad = 60 MW and P rad = 40 MW cases. When P rad = 60 MW the evolution of plasma properties is much less regular than when P rad = 40 MW. When P rad = 60 MW, excursions to high density levels at low but still outwardly directed net deuteron flux are an indication of a looming radiative collapse. As the fusion power and the power transferred from the plasma core to the SOL decrease with lower density, the plasma is more likely to become completely detached at P rad = 60 MW and runs into an uncontrollable regime that might be comparable to a density limit (dashed curves in dashed/blue colour in figure 14, on the left).
Contrarily, the cases with the lower P rad = 40 MW seem to behave much more regularly. A density limit does not appear. However, the outer target in this case does always reattach at pellet ablation (see table 1) for all feedback-controlled pedestal densities, and thus the outer target may be insufficiently protected at times of complete reattachment due to the substantial increase in heat flux after the injection of a pellet, if P rad is lowered.
The fact that n e, spx and T e, spx may follow a cyclical pattern rather than a monotonic function during a pellet injection cycle highlights the important fact that such a system can in principle be controlled by external actuators. From the controllability point of view, it seems to be advantageous to radiate less power in the SOL, as cyclically more stable conditions could then also be established at lower core plasma density and lower fusion power. It is the area enclosed by the cyclical pattern in figure 14 which has to be optimised (not necessarily minimised) to ensure a controllable scenario.
The radial particle flux Γ perp at the separatrix is usually not directly available from measurements. Instead, other parameters like neutral pressure in the divertor (from fast pressure gauge measurements at various locations) or shunt measurements of the divertor temperature [39] could be available during an ITER discharge. As the pumping capability in ITER will be restricted to ≈200 Pa m 3 s −1 or ≈ 1 · 10 23 particles per second, the estimate of the fuelling burden due to pellet injection deserves special attention. A minimal neutral flux to the pump of ≈ 0.1 · 10 23 s −1 between pellet injections is reached in the JINTRAC simulations. In the post-pellet phase, however, the limits of the pumping system might be temporarily exceeded, with an estimated maximum pumped neutral flux of ≈ 1.1 · 10 23 s −1 in the case of P rad ≈ 40 MW and ≈ 2.1 · 10 23 s −1 if P rad ≈ 60 MW. However, this result has to be interpreted with caution, as a simplified model for the geometry of the pumping structure has been used in the simulation. Other uncertainties in the model, such as transport assumptions, must also be taken into consideration. Furthermore, as already shown in [40], the fuelling efficiency and thus the fuelling burden to the pump strongly depend on the pellet plasmoid drift as well as on the effect of pellet-induced perturbations on the ELM pattern.

Summary and conclusion
An analysis of the ITER baseline scenario in the flat-top phase has been performed with the integrated modelling code suite JINTRAC, focussing on the dynamics of particle fuelling by pellet injection for the evaluation of possible ITER operational regimes and limitations that may be encountered.
From parametric sensitivity tests performed with EDGE2D-EIRENE for time-averaged steady state H-mode conditions (without ELMs and with no discrete pellet ablations), employing neon seeding to control power dissipation in the SOL, it was found that impurity seeding is vital for the mitigation of the heat flux to the target plates for an ITER H-mode scenario, assuming at least that P rad > 40 MW.
The dynamic evolution of plasma properties during a pellet injection cycle has been modelled with the JINTRAC code package. The predictions have been scanned for various levels of density using density feedback control, pellet particle E × B drift and SOL radiation P rad .
Cyclically stable plasma operational regimes could be found. For the medium density case with 50% of plasmoid E × B drift velocity and P rad = 60 MW, a quasi-stationary regime could be obtained with H 98y ≈ 0.95, P fus ≈ 450 MW, Q ≈ 13, and W th ≈ 330 MJ. For a similar case where P rad = 40 MW, SOL radiation quasi-steady state conditions were found at a lower density, with H 98y ≈ 1.0, P fus ≈ 420 MW, Q ≈ 12.5, and W th ≈ 310 MJ. However, constraints on divertor operation in terms of the maximum heat flux to the target plates and required detachment conditions for divertor heat flux control, in combination with large fluctuations at pellet ablation times, narrowed the operational space considerably. The JINTRAC integrated model predicts a rather different picture of the divertor detachment in ITER scenarios in the nuclear DT-phase. With α-particle heating from fusion processes included, an increase in density increases P fus and thus also P SOL in the relaxed state between pellet ablations. Hence, for a given P rad , the level of plasma detachment is degraded with increasing (core) density. By increasing P rad and thus power dissipation in the SOL divertor detachment can be regained, but only at the expense of a smaller operational window, i.e. the avoidance of an inverse density limit, (i.e. a looming radiative collapse) when lowering the core density through the feedback mechanism.
During pellet ablation the transient increase in P SOL is not due to the increase of the α-particle heating P fus as the αparticles thermalize on longer time-scales, but is instead due to the transient increase of convective heat flux with a local increase of density at the edge. Depending on the level of P rad , the peak target heat flux can be sustained at <10 MW m −2 . However, during individual pellet ablations the throughput capability of the ITER pumping system may be transiently pushed towards or beyond its limits by a factor of two.
It was also found that without the E × B drift of the ablated pellet material, high performance plasma operation would not be feasible with the present ITER design of the fuelling system. An improved understanding of pellet ablation is required, including pellet-ELM trigger mechanisms, to reduce uncertainties in predictions of their influence on the fuelling properties in ITER plasmas.
The simulations have revealed the importance of an integrated modelling approach with explicit consideration of transient events (avoiding time-averaging approximations) such as the injection of pellets to obtain a more realistic assessment of foreseeable scenario conditions that are in full compliance with machine and operational constraints. Because of the various simplifications that have been made in the presented simulations, the results can be considered to be valid only in a qualitative sense.
The JINTRAC integrated code is continuously being extended to improve the quantitative accuracy of the modelling predictions. The implementation of an improved discrete ELM model (including an improved inter-ELM pedestal transport model) is foreseen as one of the next steps to improve the integrated scenario modelling capabilities.
Explicit modelling of impurity transport to account for noncoronal effects should also be included as, for example, the balance of ion-thermal and friction forces may drive impurities from the divertor towards the core region, depending on the parallel T-gradient strength that depends on divertor density. In the short phase of reattachment during a pellet ablation phase, the thermal force may overcome the friction force, leading to a redistribution of impurities towards upstream. The simple mockup intoduced in the presented model to impose a constant impurity radiation loss P rad does not reflect this appropriately.
Futhermore, cross-field drift effects on plasma flows have been completely disregarded in the presented model. The inclusion of drifts in the edge/SOL plasma model usually helps in the understanding of divertor flow asymmetries in existing tokamak devices. It has been shown that drift effects seem to play a role at least in small-scale devices (refer to Alcator C-mod [41], Asdex-Upgrade and JET [42]) and it is important to understand their relevance for the validation process of state-of-the-art edge codes like SOLPS-ITER or EDGE2D-EIRENE. However, the role of drifts may be overlaid by other effects like neutral transport and it is currently being debated whether drifts play a role in large devices like ITER in (partially) detached divertor condtions. However, in view of the presented work, transients like pellet ablation and subsequent reattachment for a short period may make SOL flows more susceptible to drifts; that remains to be investigated in the future. Technical difficulty in introducing drifts in edge codes is also a common experience amongst modellers who tend to use time-steps being one order of magnitude lower to achieve code convergence. It is again a question for the pilot of the code to find the right balance between the choice of physically relevant effects in the integrated model and a reasonable run-time of the code.
The segregation effect of D and T particles and its impact on core fuelling, transport, and α-particle heating is not addressed in the current version of JINTRAC but is being investigated (see [43] for recent developments). Figures A1 and 2 show the results, assuming a higher gas flux rate of 1.4 · 10 23 s −1 in EDGE2D-EIRENE compared with a lower assumed gas flux rate for SOLPS4 of 0.5 · 10 23 s −1 . In this test, in neither EDGE2D-EIRENE nor SOLPS4 was any account taken of the molecular ion conversion (required for MAR), or neutral viscosity, or molecular elastic processes. The main difference in the model is in the divertor pump structure: unlike EGDE2D-EIRENE, SOLPS4 had two extra semi-transparent surfaces located under the D-shaped dome, in place to steer the neutral conductance between inner and outer divertor legs and affecting the total pumping efficiency. Therefore in EDGE2D-EIRENE a higher gas flux rate is necessary to achieve a similar neutral pressure in the divertor and thus to reproduce similar upstream and downstream parameters for density and temperature (see figure A1). Some discrepancies are observable for low field side target density and temperature close to the profile peaks. The heat and particle fluxes at the targets as seen from figure A2 are also in the right area (within 20-30% of their peak values) and the match for the fluxes is somewhat better than for n e and T e . We conclude that one can indeed use the EDGE2D-EIRENE model to mimic the previous SOLPS4 ITER L-mode scenarios, and that there is a way to reproduce similar divertor conditions (neglecting some details in the resulting profiles) when a higher gas flux rate is used in EDGE2D-EIRENE (~factor 3) compared to SOLPS4.
It is possible to add the extra elastic and molecular ion conversion processes in EDGE2D-EIRENE as well as a Figure A1. Electron density and temperature profiles at mid-plane and targets for the ITER L-mode scenario, comparing SOLPS4 (red) and EDGE2D-EIRENE (blue). modified divertor pump structure. This would improve the EDGE2D-EIRENE/SOLPS4 comparison by making the increase of the gas flux rate unnecessary. However, EDGE2D has a special treatment of the source terms implemented which allows the code to speed up significantly. In EDGE2D, source terms from atomic and molecular interactions with the plasma are used in a linearised form, evaluating rather neutral atomic and molecular average densities and temperatures on the EDGE2D fluid code side [17,18]. This in turn has a smoothing effect on the source terms derived by EIRENE which are generally affected by Monte Carlo noise. The consequence is a positive effect on the conv ergence rate and thus a general speeding-up of the coupled EDGE2D-EIRENE code, as is not necessary to recalculate the neutral distribution on the 2D grid at every fluid time step. EIRENE is called every ten EDGE2D fluid-time steps, which turned out to be a reasonable number in case of time steps of the order of 2 · 10 -6 s or less. The time step itself is adjusted by EDGE2D automatically, to avoid excessively strong changes in any of the source terms. A high level of Monte Carlo noise, usually seen in the ion momentum and ion energy source terms [5], is partly averaged out by the linearization procedure explained above. To ensure the particle balance between two iterative fluid time steps, the source terms are rescaled to match the integral particle source from the previous EIRENE call.
Adding extra processes like elastic or MAR on the EIRENE side is generally possible; in a dedicated code exercise it was shown that EDGE2D-EIRENE can reproduce the SOLPS4 plasma profiles (for cases where neutral viscosity can be regarded as irrelevant), assuming the same level of gas flux and a similar divertor pump geometry [44]. For this the aforementioned source linearization scheme in the coupling has to be turned off. A more detailed report is given in [23]. However, for the JINTRAC results presented, the source linearization has been kept in place to improve numerical stability and speed. For the results presented in this paper it is always assumed that the effect of extra molecular processes on the neutral pressure in the divertor can be mocked up just by increasing the gas flux, as explained above.