Simulation of neoclassical heavy impurity transport in ASDEX Upgrade with applied 3D magnetic fields using the nonlinear MHD code JOREK

Understanding and predicting impurity transport in tokamaks is necessary to ensure high plasma performance. Specfically, heavy impurities, such as tungsten (W), are very effective in radiating energy away from the plasma. Impurities can come from the limiter startup of the plasma or sputtering of the plasma-facing-components (PFCs) from which it is then transported. An accumulation of heavy impurities in the confined region of the plasma (the core) lowers the plasma temperature decreasing its performance and can lead to a disruption of the plasma.

Transport of impurities has long been studied^1,3–8 with increasingly complete theoretical and numerical models. In small to medium-sized tokamaks (MSTs), the effective (inter-ELM) radial impurity transport is typically directed inward. This then results in a peaked impurity profile, i.e., impurities concentrate in the core. When operating in H-mode, Edge Localized Modes (ELMs) typically occur. Although these instabilities are unwanted in MSTs ELMs have the favorable effect of flushing out the impurities during the ELM. Larger tokamaks, like ITER, are likely to form a hollow impurity profile without ELMs, i.e., impurities concentrate in the edge outside the confined plasma. The reason for this is the difference in the collisional neoclassical transport. Bigger tokamaks have a steeper temperature gradient in the pedestal, while density gradients are shallower since a high separatrix density is required for divertor power exhaust. This results in stronger temperature screening for the impurities. This shifts the balance from effective inward convection to effective outward convection. It has been shown that for hollow impurity profiles ELMs will suck in the impurities into the confined plasma core.¹ The combination of the heat load impacts of bigger type-I ELMs and the inward transport even during smaller ELMs drives a need for ELM-controlled or ELM-suppressed H-mode plasmas.

It is possible to mitigate and suppress ELMs by applying external 3D magnetic fields, Resonant Magnetic Perturbations (RMPs). The RMPs cause ergodization of the edge magnetic field and 3D-structures instead of the axisymmetric flux surfaces. Due to the ergodized edge, the plasma pedestal decreases, and heat and particle fluxes increase and the plasma becomes inherently 3D. In experiments in AUG, an enhanced flushing out of heavy impurities is seen due to the use of RMPs for ELM-suppression. However, there is no analytical model that predicts the effective radial impurity transport. Understanding and correctly predicting the impurity transport can have a major impact on how ITER (or other future tokamaks) can operate.

There is increasing effort in numerically studying W transport, either during ELMs^1,9 or in RMP ELM-suppressed plasmas.^10–12 The RMP ELM-suppressed studies all show enhanced transport on a static background. However, there are still limitations to the models for understanding the impurity transport through the pedestal.

JOREK is a fully implicit 3D non-linear extended MHD code for realistic tokamak X-point plasmas.^13–15 The JOREK code has typically been developed and used for studying MHD instabilities such as ELMs and disruptions.^16–18 The present version of JOREK includes some divertor physics in terms of sheath boundary conditions (Mach-one velocity and sheath transmission factors) and a simple neutral fluid model (reflection of outgoing ions as neutrals, ionization, recombination, and radiation). This model has been extensively applied for simulations of ELMs¹⁸ and ELM control by RMPs¹⁹ and by pellet triggering.²⁰ JOREK has recently been extended with a kinetic neutral model²¹ to improve upon the description of the Scrape-Off-Layer and divertor.

The aim of this paper is to describe benchmark and apply an improved physics model to investigate and understand W transport in RMP ELM-suppressed pedestals. As W transport is determined by a multitude of factors—not all of them understood yet with applied 3D fields—as basis, we select one scenario from ASDEX Upgrade (AUG) and form an initial understanding of the effect of RMPs on W transport.

This paper is structured as follows: we first summarize the JOREK models used. Then, in Sec. III, we present the applicable neoclassical theory and validate our model with said theory. Section IV describes the simulation set-up and generally how RMPs affect the plasma. In Sec. V, a scan over RMP current is performed and its relation to effective W diffusion is discussed. Section VI shows the effect that the electric fields and collisions have on radial transport in one scenario. Section VII presents the change in permeability of the pedestal between an axisymmetric plasma and ELM-suppressed plasma with an identical pedestal. It shows enhanced transport from both inside and outside the separatrix. In Sec. VIII, the full 3D picture is explained, where the possibility of W trapping in potential wells is shown and discussed. Section IX introduces experiment-like tungsten laser-blow-off sources to simulate transport from the wall to the core. Due to a lack of suitable diagnostics for tungsten in the pedestal, this section also investigates the use of argon as a substitute. Finally, Secs. X and XI are the discussion and conclusion.

This section briefly summarizes the used MHD model in JOREK in the reduced MHD form to provide a background for this work. A more complete explanation of the JOREK model can be found in Refs. 13–15. The development of the particle framework and coupling projections of moments of the particle distributions to the finite elements can be found in Refs. 1, 21, and 22. The development and recent work on Free-boundary ELM-suppression simulations can be found in Refs. 23 and 24.

Here, ψ is the poloidal magnetic flux. The velocity stream function u is defined as $Φ / F 0$⁠, the electric potential $Φ$ divided by F₀. The reductions give us seven equations with the seven variables: poloidal magnetic flux ψ, velocity stream function u, toroidal current density j, toroidal vorticity w, mass density $ρ = m i n i + m e n e$⁠, MHD temperature $T = T i + T e$ and parallel velocity $v ∥$⁠. There is also the possibility to run with a 2-temperature model in JOREK, but is not used in this work.

For the initial equilibrium, JOREK has an internal 2D Grad–Shafranov solver, which uses the same finite element representation as for the further time-dependent simulations. This guarantees the initial state in the MHD equations satisfies the initial equilibrium force balance on the discrete level. The solver utilizes pressure and FF $′$ profiles as a function of the normalized poloidal flux, and a list of (⁠ $R , Z , ψ$⁠) points defining the boundary of the computational domain. Keeping ψ_bnd constant over time means fixed-boundary MHD.

To model the interactions between external conducting structures (e.g., in-vessel coils, resistive wall) and the plasma, JOREK is coupled to the vacuum-field code STARWALL.²³ STARWALL couples to JOREK via boundary conditions of the plasma current density and poloidal magnetic flux. This replaces the otherwise applied Dirichlet boundary condition (corresponding to an ideally conducting wall) by natural boundary conditions. We denote this “free boundary treatment.” Note that for this work, we have chosen to keep the axisymmetric component to be fixed to avoid dealing with vertical plasma instabilities or position control. However, we treat the n > 0 modes with this free boundary approach such that the 3D perturbation fields on the computational boundary are self-consistent.

On top of the non-linear (fluid) MHD model, JOREK is extended with a Monte Carlo particle tracker, which can be two-way coupled to the MHD fields. The model utilizes a full-orbit pusher, including ionization recombination, effective lines, and continuum radiation based on the background fields. The effective collisional-radiative rates come from the OpenADAS database.²⁵ No assumptions are made on the impurity charge distributions as all charge states are modeled. The evolution of the charge state of each marker particle is governed by effective rate coefficients of ionization and recombination.

Once initialized, marker particles stay alive between each MHD fluid time step, until they are removed by a reaction or by hitting the wall. This ensures a continuous and consistent particle distribution over time.

Arbitrary moments of the particle distribution function can be projected onto the JOREK finite-element representation, either for diagnostic purpose or for feedback onto the plasma.

This section presents the collision operator for binary collisions, as implemented in the JOREK code. Then, per subsection, a component of collisional neoclassical transport is explained and a benchmark of the collision operator is performed.

As we intend to study impurity transport in the H-mode pedestal with applied 3D-magnetic fields, typical flux-surface-average neoclassical transport is not suitable. Here, the combination of steep gradients, $ρ z ∼ L t = T | ∇ T | , L n = n | ∇ n |$⁠, ergodic magnetic field lines, and the choice that JOREK utilizes fully implicit time stepping indicate a fluid neoclassical model is not the right approach.

We employ the hybrid fluid-kinetic model previously developed in JOREK.^1,21,22 Here, impurities are modeled as full-orbit marker particles representing many physical particles and are pushed based on the local MHD-fields with a Boris integrator.^26,27 In this approach, neoclassical transport should arise as a collective effect due to many binary collisions between the impurity marker particles and the background plasma. Therefore, we need to understand the basics of neoclassical transport and the adoption of a collision operator.

Collisional Neoclassical impurity transport arises due to the magnetic topology of the plasma and is driven by drift-orbits, frictional and viscous forces. For now, only collision with background main ions will be considered. The impurity ions undergo (small angle) Coulomb collisions. Viscous forces are ignored as sufficiently high background collisionality is assumed. In a tokamak, the magnetic topology causes enhanced diffusion. A perpendicular density gradient causes inward convection, while a perpendicular temperature gradient results in outward convection.

To verify and validate the collision operator several simulations have been performed. In  Appendix A, the test prescribed in Ref. 30 is presented and shows good agreement with the theory. These tests show the JOREK particles describe well the thermal forces on the particles and can simulate the temperature screening effect.

Second, in this section, a neoclassical diffusion benchmark is presented. Here JOREK test particles are simulated on a background circular deuterium plasma with characteristic parameters shown in Table I. The plasma is chosen such that the high-collisionality neoclassical transport model is applicable.

The W ions are in the Pfirsch–Schlüter regime, $ε 3 / 2 ν z * ≥ 1$⁠, and the inverse aspect ratio ε is low enough that inner-outer asymmetry should be negligible for the diffusion.³¹ 

The W ions are initialized as a spatially Gaussian profile at normalized poloidal flux values listed in Table I with a half-width half-maximum of $0.07 ψ n$⁠, and with the same temperature as the plasma background. The particle pusher (Boris pusher) time step, $Δ t kinetic$⁠, and collision time step, $δ t coll$⁠, are chosen to resolve the particle orbits and the neoclassical collisions.

Using the guidelines of Homma, the advised gyro-time step and binary collision model (BCM)-time step are calculated as a function of the background plasma. Figure 1 shows the advised time steps for resolving the gyro-orbit t_gyro and the binary collisions t_BCM as a function of the plasma temperature for the range of all charge states for W ions. The background parameters are the same as used in the validation benchmark as shown in Table I. The red dots correspond with T = 250 eV with Z = 19 the charge state of W. This is the most probable charge state at that temperature in coronal equilibrium. t_coronal shows the effective time step needed, where W is in coronal equilibrium with the background plasma. The number of binary collisions needed per particle time step $n coll = t gyro t BCM$ if the maximum allowable t_gyro is used. The right axis in Fig. 1 shows the normalized collisionality of the electrons $ν e * ε 3 / 2$ and W ions $ν z * ε 3 / 2$⁠. $ν z * ε 3 / 2$ is above the PS-limit in all the locations used in this benchmark.

Note that from Fig. 1, it is relatively simple to estimate a proper time step for different scenarios. $t BCM ∝ n e − 1$ and $t gyro ∝ B − 1$⁠. In typical tokamak regimes, the coronal equilibrium is practically independent of n_e in steady-state.

In this subsection, we investigate the convergence of the diffusion coefficient as a function of the number of collisions per time step, and verify the collision operator to neoclassical diffusion.

Figure 2 shows the convergence of the Diffusion coefficient as a function of the collision frequency. The dotted vertical line represents the advised maximal collision time (or minimum collision amount) in order to have a converged solution. Figure 2 shows that indeed the Diffusion coefficient decreases below this threshold. However, the error is relatively small, even when $( t gyro t BCM )$ is 2 instead of the advised 7 in Fig. 2. Oversampling the collisions does not increase diffusivity.

The advised $δ t coll$ is used for the neoclassical diffusion benchmark. At each initial position, the particles are simulated for $≈ 50$ ms. It takes 5–10 ms before the expected diffusion behavior is established. From the profile widths over time from 15 to 45 ms a fit is used to determine the diffusion coefficient as shown in  Appendix B. The resulting normalized diffusion coefficient $D / D C L$ as a function of the safety factor q is shown in Fig. 3. The horizontal error bars on the JOREK diffusion coefficient indicate the FWHM (Full-Width Half-Maximum) of the initial particle distribution. The solid lines indicate the classical (CL), Pfirsch–Schlüter (PS), and summed CL + PS theoretical diffusion coefficient. The dashed line is a fit of the JOREK results.

The JOREK results agree well with the predicted $D neo = D C L + D P S$ as in Eq. (6). The fit shows JOREK has on average a 4% higher diffusion coefficient. It is suspected that the slight increase in D is due to the finite width of the JOREK particle profiles. This means part of the particles sit at higher q values and diffuse faster. This is also the case for the neoclassical transport tests in Ref. 35, where the neoclassical transport parameters are skewed to slightly higher q values.

Sections III A and III A 1 present a benchmark on a circular plasma with a flat density and temperature profile. Neoclassical collisions then only result in extra diffusion. This subsection presents the same circular plasma but with the effect of a density gradient and the effect of a temperature gradient. The former causes a neoclassical inward pinch of the impurities, the latter a neoclassical temperature screening.⁴ 

The results of these tests are shown in Fig. 4. It shows that the obtained average radial velocity as a function of the impurity charge states and a theoretical estimation for three different locations in the plasma. The solid black line indicates the theoretical radial velocity at the mean of the initial particle distribution. The other dotted and dashed lines are the theoretical predictions on the half-maximum of the distributions. “Half min” is at the half-maximum of the distributing for the minimum radius. “Half max” is for the half-maximum at the maximum radius.

Figure 4 shows that the JOREK results are satisfactorily close to the neoclassical theory with approximately linear scaling with the impurity charge state Z_z. On the low Z side aligns better with the half-max line where q_safe is higher. Z_z = 2 results in a higher IWP velocity than expected. Initially, the Z = 2 transport is slower but on average it results in the same transport as for Z = 3. Details of the average position of time are shown in  Appendix C.

For the temperature screening effect, a series of tests as proposed in Ref. 30 have been performed and shown in  Appendix A. These tests show that the particle model can correctly simulate the diamagnetic thermal force and the temperature screening effect in a cylindrical plasma. For toroidal geometry, again we do the same as in Fig. 4 but with a fixed temperature gradient. Figure 5 shows the result of the TSE radial velocity as a function of the impurity charge state. As with the IWP, a linear scaling with Z is observed for the TSE.

In the previous comparisons—Figs. 4 and 5—fixed convective coefficients were used; K = 1 and $H = − 1 / 2$⁠. The ratio of these, $H TSC = H / K$⁠, is the effective temperature screening coefficient. From the NEO code³⁸ and several models,^3,4,37 we know there is a dependency on the collisionality, the mass, charge state, and pressure anisotropy. The simplest model in the Pfirsch–Schlüter model in the small ε limit for $Z z ≫ Z i$ (small ε PS-model) where there is no dependency and $H TSC = − 1 / 2$⁠. The Hirschman–Sigmar–Wenzel model includes a numerically calculated dependency on normalized collisionality $g = ν i * ε 3 / 2$⁠. In a poloidally symmetric, homogeneous, non-rotating, trace impurity (⁠ $α ≡ n z Z z 2 n i Z i 2 → 0$⁠) limit, the HSW-model reduces to $H TSC ≈ − 1 2 + 0.29 0.59 + 1.34 g − 2 − Z i Z z$⁠. At $g ≫ 1$⁠, ion-electron collisions become non-negligible³⁷ and at low collisionalities, there is a reduction of the thermal force due to Banana-Plateau (BP) contributions;³⁷ these result in stronger dependencies on g. FACIT,³⁷ an analytical model that can reproduce transport coefficients from the NEO code, is here assumed to be the “ground truth” for collisional neoclassical transport.

To better understand what the effective temperature screening coefficient (TSC) is, and how it compares with theory, a scan of the average radial transport has been done as a function of the ion collisionality. Figure 6 shows the H_TSC for $W 16 +$ (Z = 74, A = 183.84) over a range of ion background collisionalities at constant pressure. The HSW-model is shown for $W 16 +$⁠. To obtain a fuller indication of H_TSC, FACIT results are shown with the average charge state of the coronal equilibrium at each collisionality.

It shows that our used collision model from Homma does not have a (clear) dependency on collisionality. Herein, it aligns well with the large aspect ratio (small ε) approximation Pfirsch–Schlüter model but with the HSW-model correction of the finite impurity charge state. Both at the very high and low collisionality, reductions of $| H TSC |$ are not reproduced.

From Fig. 6, we conclude the convectional Homma collision operator, with Spitzer–Härm heat flux,³⁵ operates well in the complete pedestal range in AUG. However, it does not include dependencies on collisionality both on the low and high collisionality range. For this purpose, this collision operator is sufficient for the scope of this research. With neoclassical diffusion, neoclassical inward pinch (IWP), and temperature screening effect (TSE), our improved JOREK particle model can simulate neoclassical impurity transport up to satisfactory accuracy.

The top of the ITER pedestal will have lower collisionality than AUG, implying our adopted collision operator overestimates the temperature screening at the top of the pedestal. The implications of using this collision operator for ITER-like plasmas are further discussed in Sec. X A.

In the previous sections, the particle model has been successfully bench-marked against neoclassical theory. The limits of the models are known, and therefore we'll focus on transport in the pedestal. This section presents the plasma scenario and simulation set-up for the ASDEX Upgrade simulation with applied 3D fields. The simulation is modeled after shot numbers as shown in Table II.

For ELM suppression in AUG, typically the current in the RMP coils I_rmp is in the range of $0.9 − 1.2$ kA and the other plasma parameters in Table II remain similar in other RMP experiments in AUG. All listed shots in Table II have toroidal n = 2 periodicity for the RMP perturbation. There is a $Δ ϕ = 90 °$ phase shift between the upper and lower row of coils.

The base scenario for the JOREK plasma is the same for all the cases, but is modeled after #34548. Fits on experimental diffusivity profiles $D ( ψ N )$ and heat conductivity $κ ( ψ N )$ together with an estimated auxiliary heating profile and particle source are used to maintain the experimental density and temperature profiles constant in time. A toroidal momentum source is added to match the experimental toroidal rotation of the plasma.

The grid extends from the plasma core almost to the first wall. For simplicity, a flux-aligned wall is chosen which only deviates a couple of centimeters from the real wall around the mid-plane. A simplified geometry for the divertor is used. Without neutrals, and trace impurity particles, this does not change the upstream plasma condition.

In Secs. V–IX aspects of impurity transport in the simulations are addressed. Although there is some global knowledge and some dedicated experiments focused on W transport with applied 3D fields for ELM-suppression, there are no direct measurements of W in the pedestal. Ergo, a one-to-one comparison with the experiment is not (yet) feasible.

As impurity transport through the pedestal is sensitive to many parameters, a generalization of impurity transport with applied 3D fields is complicated and not the goal of this study. We, therefore, focus on some general aspects that will influence cross field transport. For the chosen scenario, a scan in RMP current is performed. The effect of the Electric field and collisions on radial convection is studied. Enhanced permeability of the pedestal due to ergodization is observed. Further, it is found that tungsten could be trapped in electric-potential wells outside of the separatrix.

Finally, we present a tungsten laser-blow-off simulation to involve more accurate tungsten sources in the simulation.

Although, in our study, the focus lies on tungsten transport, it is helpful to discuss some general effects of RMPs on the plasma in the simulation and experiment to understand how these can affect tungsten transport.

The in-vessel coils (RMP coils) generate a small non-axisymmetric perturbation compared to the toroidal field (⁠ $δ B / B toroidal ∼ 10 − 3$⁠).³⁹ The plasma responds to these perturbations and has a particularly large effect on resonant flux surfaces, where they breakup the field lines and induce magnetic islands.²⁴ Due to the RMP effects, the pedestal height decreases, especially the density pedestal is decreased due to density-pump-out effects. There is still ongoing research to obtain a complete view of the mechanisms needed for achieving ELM suppression due to resonant magnetic perturbations. There are several hypotheses: One important ingredient for the ELM suppression is the changes in pedestal structure due to the RMP-driven transport.⁴⁰ In ASDEX Upgrade, there are observations of broadband turbulent fluctuations and the existence of a quasi-coherent mode which are thought to reduce the pressure gradients in the pedestal.⁴¹ Gyrokinetic turbulence simulation of the H-mode pedestal in the presence of RMPs shows the erosion of zonal flows and an increase in the electron-scale transport.⁴² Another important ingredient is toroidal mode coupling. Above a certain threshold, the toroidal coupling of the unstable edge modes with the static mode induced by RMPs is large enough to lead to saturation.^43,44 Thus, higher-order toroidal modes in the edge, responsible for the ELM-crash, cannot unstably grow beyond the saturated state. The electric field becomes 3D due to the RMP as well and is a part of the self-consistent 3D MHD solution.

The effects of the RMPs on the plasma are visualized in Fig. 7. It shows a poloidal cross section of the electron density where perturbations in the edge become visible. The Poincaré plot (white dots and lines) on top of the density show the magnetic field structure. The toroidal part indicates the magnetic poloidal flux perturbation. It is important to understand the extent to which magnetic perturbations penetrate into the plasma and until where the effects of RMPs occur. Despite the finite resistivity—with which the plasma would not screen the perturbation on longer time scales—plasma rotation can induce mirror currents on resonant flux surfaces, which counteract the local resonant component of the perturbation and prevent the penetration of magnetic islands.^45,46

Furthermore, due to the perturbations, instead of two axisymmetric plasma legs forming, several lobes are induced by the destruction of the separatrix. This causes multiple non-axisymmetric strike lines on the divertor.

This section. summarizes the particle physics included and how the simulations are set-up for Secs. V–IX.

Tungsten is treated as a full-orbit trace particle. The Boris pusher is used to update the particle location and velocity; here the local magnetic field and electric field vectors are included and no assumption is made on the existence of a flux surface. It thus includes, gyration, and grad/curvature-B drifts due to the magnetic field; and $E × B$ drifts and parallel E acceleration due to the electric field. The charge state of the marker particle is governed by a Monte Carlo approach with the ionization and recombination probabilities from effective rate coefficients. Effective radiation is included, but does not affect the plasma. The Homma collision operator is used with neoclassical friction [as presented in Eqs. (2) and (3)], i.e., parallel flow and temperature gradients shift and distort the local Maxwellian from which is sampled. Here, only impurity-ion collisions are taken into account.

From the same scenario, four free-boundary simulations are set-up with different RMP currents (⁠ $I RMP = [ 0 , 0.5 , 1 , 1.22 ]$ kA). The perpendicular transport fitted for the I_RMP = 1 kA case remains unchanged for all the RMP cases. In Sec. V tungsten is initialized in location 1, as shown in Fig. 8, for all four RMP currents. In this set-up, we probe the tungsten diffusivity in the pedestal.

To understand the role of the E-field and collisions in Sec. VI, tungsten is also initialized in location 1. Three identical simulations are performed with exception of including/excluding the following physics in the particle pusher: either with electric field plus collisions; without electric field plus with collisions and with electric field without collisions.

For understanding enhancement of pedestal permeability in Sec. VII—both transport from outside to in, and vice versa—, tungsten is both initialized in location 1 and 2 in Fig. 8. This is a simplified way to understand transport through the pedestal without incorporating a global tungsten profile with realistic sources and sinks. In each location, two different background plasmas are utilized to compare the tungsten transport. To solely investigate the implications of the 3Dness and ergodization of the plasma, transport coefficients are altered ad hoc when the RMPs are switched off such that, the density and temperature pedestal remain identical. The fields, currents, and flow do evolve self-consistently. Then, tungsten trace particles are simulated on top of the axi- and non-axisymmetric background.

Finally, a more experimentally relevant tungsten source is added by means W-LBO (tungsten Laser-Blow-Off) in Sec. IX. Here, a relatively cold W distribution is injected from the simulation edge at a poloidally and toroidally localized point.

A scan over a range of RMP currents, both below and above the ELM-suppression threshold, has been performed. This includes two sub-threshold cases, I_RMP = 0 and 0.5 kA, where I_RMP = 0 is the normal axisymmetric case; and two ELM-suppressed cases with I_RMP = 1.00 and 1.22 kA. The plasma responds to the external fields in a multitude of ways (e.g., kink response, ergodization, and shielding). Each response might contribute to altering the impurity transport across the pedestal. In Fig. 9 we present an effective tungsten diffusivity, $D W , eff$⁠, in the pedestal region as a function of the RMP currents.

Figure 9 (main plot) shows a non-trivial increase in $D W , eff$ as a function of the RMP current (red stars). The solid blue line indicates what change in $D W , eff$ is expected from the changed n, T-profiles due to the RMP current (i.e., neoclassical diffusivity). Although, especially the density at the pedestal top decreases due to RMP density-pump-out, it is not expected that the diffusivity changes more than 10%. However the effective diffusivity increases by about a factor of 2 going into the ELM-suppressed regime, which cannot be explained due to the changing background pressure (Notethat the vertical axis starts at 2.5).

The insets in Fig. 9 show a corresponding Poincaré plot of the magnetic field lines for each point in RMP current in the normalized flux range $0.93 ≤ ψ N ≤ 1$⁠. Note that there is also a kink response in the toroidal direction which is not captured in the inset plots. The black dashed–dotted lines in Fig. 9 (insets) are where the W particle distributions are initialized, with a standard deviation of 0.012ψ_N. Although at 0.5 kA, there is a kink response and some ergodization at the edge, the diffusivity only increases marginally. Once in the ELM-suppressed regime, where the particles see a strong kink, magnetic islands and stochastic field lines, it significantly increases the diffusivity.

Note here, that although the magnetic field structure at I_RMP = 0.5 in Fig. 9 (insets) looks strongly ergodized, the field line connection length is still long (i.e., one to tens of kilometers). Collisionless diffusion due to the magnetic field line ergodization is only a small effect, as is discussed in Sec. VI. Furthermore, the background ion diffusivity does not change throughout the different simulations. The perpendicular diffusion and heat conduction coefficient have been fitted such that the experimental profiles for the density and temperature are obtained.

Beyond effective diffusion, it is important to understand the role of various physical components on the radial (i.e., cross field) convection. In this section, we investigate the roles of collisions, the electric field, and their synergy. These simulations are run with I_RMP = 1 kA and thus in the ELM-suppressed regime. Only physics mechanisms in the particle framework are altered, the plasma background remains unchanged.

Figure 10 illustrates the effects of adding either collisions, or the electric field in the particle pusher, or both. The figure is like a “streak camera” with time on the vertical axis and distance from the separatrix (projected on the OMP) and (⁠ $Δ r separatrix$ at OMP), on the horizontal axis. The color indicates the flux surface averaged component of the normalized W density. Figure 10 (top) shows the density profiles of a time slice at 24 ms. Figure 10(a) is the full simulation including the effects of electric fields and collisions on the particles, (b) with collisions but without the E-field effects, and (c) without collisions, but with the E-field.

Without collisions, Fig. 10(c), the transport along ergodic field lines is comparatively slow. Even, the E × B drifts do not add significant transport here. Figure 10(b), including collisions but w/o E-field, we obtain collisional neoclassical transport, but on ergodic field lines. It shows that convection drives the distribution both in- and outward. The two branches find a new region with neoclassical flow cancellation. These branches are located respectively at the bottom of the pedestal or in the SOL, and at the top of the pedestal. Although, the convection speed inwards appears slower.

Figure 10(a), including both the effects of the E-field and collisions, we observe the importance of including the effects simultaneously. The collisional neoclassical transport drives the two branches, but the effective stagnation locations shift due to the electric field. The population fraction transporting up the pedestal increases by including the E-field.

Note that the results presented here are the flux-surface averaged components projected onto the OMP around the separatrix. Poloidal and toroidal asymmetries disappear in this projection. There is further no physical significance to projecting the tungsten density onto the OMP.

Although direct projection of the perpendicular-particle-flux moment of the distribution function (i.e., $∫ v → f · ∇ ψ / | ∇ ψ | d V$⁠) fluctuates strongly per time step, over longer time (i.e., millisecond) a local convection velocity becomes clear from Fig. 10. The approximate radial velocity, v_r, of Fig. 10(a) in $0 ≤ t ≤ 7$ ms is −3.7 m/s for the portion branching inwards.

Let us further inspect how the electric potential and thus the E × B profiles change in response to the RMP fields. Figure 11 shows the electric potential contours of the top part of the plasma. The E × B flow direction lies on the potential contour lines and are indicated with red arrows. The pink contour is the unperturbed (i.e., without RMPs) separatrix. The white iso-contours are the unperturbed electric potential, and the black iso-contours are the potential perturbed by the RMP fields. In Fig. 11, there are only vortices turning in one direction. Their counterparts, i.e., the vortices turning in the other direction that are expected in between the vortices shown, are compensated by the n = 0 flow pattern. Note that the density flows in and out according to the E × B flow direction.

The E × B vortices lie mostly outside of the separatrix, indicating that in Fig. 10 the E-field does not have a dominant effect. However, for W coming from outside of the separatrix, E × B drifts are probably a dominant transport mechanism in the region of the E × B vortices.

Apart from the previously described changes in diffusive transport (Sec. V) and the role of the electric field (Sec. VI), the question on how neoclassical cross field convection is altered due to RMPs is discussed in this section. To solely investigate the implications of the 3Dness and ergodization of the plasma, the background fluid transport coefficients are altered ad hoc such that the n = 0 component of the density and temperature profile remains unchanged. Flows and currents are evolved to obtain a self-consistent equilibrium. Here, we present two cases on W transport across the pedestal, one initialized outside, and one initialized inside of the separatrix. It further corresponds respectively to the bottom of the pedestal, and the middle of the pedestal.

Figure 12 shows W transport through the pedestal from outside the separatrix, both with RMPs on, and in the axisymmetric case without RMPs. In both cases, collisions and effects of the E-field in the transport are turned on. As in Sec. VI, the results are the axisymmetric component projected onto the OMP.

In Fig. 12(b), without RMPs, W on the open field lines quickly flow out to the divertor. Hence, the effective inward transport is very low and W is well-screened from the pedestal. Note, without neutrals in the simulation, the plasma operates in the sheath-limited regime which is characterized by parallel convective heat transport and a low parallel temperature gradient in the Scrape-off-Layer (SOL). In these conditions, the thermal force driving impurities upstream is low. With the RMPs, Fig. 12(a), there is significant W transport across the separatrix but the inwards transport stagnates at the bottom of the pedestal and does not lead to W transport to the pedestal top. This fast transport from just outside to just inside the separatrix is likely enhanced by the formation of E × B vortices as shown in Fig. 11.

Although, the change in permeability of the pedestal looks drastic from the figure, the radial stagnation point is located at the bottom of the pedestal. There is no strong convective transport through the pedestal into the core on this timescale.

When W is initialized in the pedestal, Fig. 13, enhanced pedestal permeability is also observed. Figure 13(a) is the same as Fig. 10(a); it shows that the branching to the top and bottom of the pedestal is much faster than without the RMPs as seen in Fig. 13(b). In the axisymmetric case, there is transport to the separatrix, but the bulk diffuses and convects significantly slower.

The two case studies (i.e., initializing W just inside, or just outside the separatrix) show that RMPs enhance transport through the pedestal. As the n-,T-profiles were kept identical, the flux surface averaged collisional neoclassical transport should be the same according to the model [Eq. (7)]. Of course, the real transport changes, the plasma becomes non-axisymmetric and there is also an electric field response to the RMPs. The combination of the non-axisymmetry and ergodization cause an effective enhanced permeability of the pedestal. The electric field response is clearly important for the transport, but without the effect of the E-field on the particles, the enhanced permeability is still observed as shown in Fig. 10.

Note that the actual effective transport and screening properties without RMPs are different than in these case studies. After turning off the RMP currents, the plasma will not only become axisymmetric again; the pedestal will also restore. Meaning the temperature and density gradients change, altering the collisional neoclassical transport. Notably, the density pedestal will be significantly steeper which causes more inward transport.

So far, we have mainly presented flux-surface-averaged quantities of transport, while in reality, the transport is strongly non-axisymmetric. Further, the previous projections were all projected onto the OMP; however, there are also strong poloidal asymmetries; strongly affecting transport as well.³¹ This section presents the 3D and non-axisymmetric transport of tungsten.

Figure 14 shows quasi-steady-state poloidal W density cross sections and contours of W density with I_RMP = 1 kA. Here, the pink contour indicates the highest density of W, gray is medium W density. The two contours bear no physical meaning besides their purpose for visualizing important aspects. W was initialized in the pedestal (location 1). A small part convects into the core which most convects to the bottom of the pedestal as in Fig. 13. Strong localized poloidal asymmetries develop, especially on the high-field-side (HFS). The asymmetries are of helical nature due to RMPs. From the gray contours it becomes apparent that the HFS poloidal asymmetries are continuously connected.

On the upper region of the plasma, the pink contour reveals a continuous toroidal ring, showing that W is not trapped in this region. However, toroidally and poloidally localized blobs form in other parts of the plasma shown in pink; notably on the low-field-side (LFS), outside the separatrix, close to the X-point. In these regions, W gets trapped in the blobs; not in a collisionless banana-orbit way, but collisionally in an electric potential well.

The 3D potential well structure becomes visible by projecting the helical potential perturbation onto a toroidal line, or onto a line perpendicular to the magnetic field. Then $δ u → · ∇ ψ ≠ 0$ and $δ u → · ∇ ϕ ≠ 0$⁠. Along such a line, the perturbation appears as a potential well. Figure 15 shows the electric potential and the W density on a toroidal line that intersects the W blob. It shows W only gathers in the well. Although it has a toroidal periodicity of 2, up to toroidal n = 16 harmonics are needed to resolve the potential structure sufficiently. A poloidal cross section also shows a potential well.

W convects into the well and, on longer time-scales, it diffuses out of the well. The occupation in the well in the experiment would depend on the W sources. In addition, there are a couple of mechanisms that might affect this quasi-static solution. First, radiation feedback on the plasma lowers the temperature and might lower the potential well increasing the W trapping. Second, at sufficient W density, self-diffusion (i.e., W–W collisions) is not included in the current model. This would limit the W accumulation. Finally, electrostatic turbulence, which is not included in our model, and can influence the extent to which such a quasi-static solution is possible.

In the ELM-suppressed regime, the electric potential forms such structures and thus will tend to accumulate W in such 3D potential wells.

Experimentally, a controlled W source can be produced by means of laser-blow-off to investigate the W transport. A laser burst evaporates a known quantity of W on a glass plate situated on the LFS around 30 cm below the OMP in AUG.⁴⁷ Although intrinsic impurity content often cannot be ignored in such experiments, radiation due to the sudden W influx provides ample signal-to-noise ratio.

In this section, a W-LBO source is added in the I_RMP = 1 kA scenario. This could give the possibility to compare to the experiment. The simulation is of very similar nature as in the previous sections. However, W is initialized at the edge of the simulation domain as neutrals with a relatively low velocity sampled from a cosine distribution. This source is poloidally and toroidally localized (⁠ $50 × 50 mm 2$⁠).⁴⁸ Due to the high collisionality of W, the initial low-temperature W is quickly equilibrated with the background plasma temperature.

Figure 16 shows the initial phase of the W-LBO simulation, where W is dragged up to the background velocity due to friction and transports to a quasi-stationary convection-stagnated state. Here, only the n = 0 projection is shown (i.e., the axisymmetric component) with a logarithmic color bar.

W does not penetrate far into the plasma before being strongly dragged along with the plasma. The strong parallel friction causes inward convection of W. In around 7 ms perpendicular convection stagnates. Again, W gets trapped into the potential well. In the first 7 ms, 65% of W is lost to the wall. After the initial convective phase, a diffusion-limited phase occurs where W is located mostly at the bottom of the pedestal and transport is slower. In about 70 ms, 95% of W has been lost.

Although in Fig. 16 it seems as if W is transported poloidally, Fig. 17 shows that this transport is mostly parallel to the magnetic field, as expected from friction with the background plasma. Note that the contours show the highest densities of W. W is not solely transported with this bulk. This can be observed in the (back) cross section in Fig. 17 and in Fig. 16.

Due to the effect of the RMPs, there are strong flow shears in the edge, thus W can locally undergo a strong increase in friction. Note that e.g., I_RMP, RMP shielding and toroidal rotation can influence this transport, as it locally will alter the background velocity and mode structure.

Although, with these simulations, transport through the pedestal can be investigated, it remains difficult to validate with the experiment. The spectroscopic tools used in AUG only measure transitions at high ionization states of W. Thus, below the ∼keV temperature range, W transport is not observed.

Due to the lack of suitable diagnostics for W in the pedestal, our model cannot be directly validated by measuring W transport in the pedestal, a set of similar simulations is performed with argon as an impurity species by introducing a short gas burst from the edge and initializing Argon in the same way as was done in Sec. VII. Argon is still a relatively heavy impurity compared to deuterium and radiates at lower temperatures than W but has the advantage that it would be measurable in the pedestal.

The simulation set-ups are the same as before. However, argon is now the impurity including all the effective rate coefficients and ionization states. This subsection investigates how similarly Ar behaves compared to W in the presence of RMPs and whether it could be used as a proxy to validate our simulations with a possible experiment.

The first simulation introduced a short gas burst similar to W-LBO. In this simulation, we probe whether injected argon in an experiment would get trapped in potential wells.

Figure 18 shows part of the Ar distribution sufficient time after the Ar-puff. Although the background is identical—and thus 3D potential wells are present—Ar does not get trapped. It also has a broader footprint in the SOL and pedestal indicating enhanced transport through the pedestal compared to W. Although Ar does not get trapped, which indicates some difference in transport, this does not yet mean argon cannot be used as a substitute.

To further understand the difference between W and Ar transport, simulations are performed where Ar is initialized at the top and bottom of the pedestal identical as in Sec. VII. Two cases on Ar transport across the pedestal, one initialized outside, and one initialized inside of the separatrix. It further corresponds, respectively, to the bottom of the pedestal, and the middle of the pedestal. It includes all effects in the model. As Ar has different ionization potentials, ionization and recombination rate coefficients, the coronal equilibrium is different for the same background plasma (i.e., different charge state probabilities).

Figure 19 shows Ar transport through the pedestal, both from the inside (a) and outside (b) of the separatrix. Both cases include the RMPs. The results are the axisymmetric component projected onto the OMP. The stagnation of radial convection just inside the separatrix occurs at the same location as for W [see Figs. 10, 12, 13(a), 19(a), and 19(b)]. This means the transport direction changes between inwards and outwards at the same radial position for Ar and W. Transport from both initial locations behaves similarly to W, but it has some important differences. Ar is more diffusive than W. Figure 19(b) shows more diffusive penetration into the pedestal. Figure 19(a) shows effectively more inward convection. Both enhanced diffusion and cross field transport for Ar is consistent with neoclassical theory.^49–51 Less Ar is driven to the flow-cancellation point just inside of the separatrix, and more is driven further inward. W encounters another stagnation point at $Δ r separatrix = − 4.2$ cm, while after 28 ms the innermost Ar peak lies at $Δ r separatrix = − 5.9$ cm and is still slowly drifting inwards with approximately 0.5 m/s. The fastest inward convection of Argon $v r , A r ≈ − 5.2$ m/s which is 1.5 m/s faster than for W.

This indicates that with applied 3D fields, Ar behaves similarly in terms of transport through the pedestal. Although the pedestal is more permeable to Ar and inward transport is enhanced, Ar could well be used to study W-like transport through the pedestal; especially because the stagnation (or flow-cancellation) point lies at the same radial position. Further, it could be used to validate our models with the experiment.

After verifying the collision operator we investigated tungsten transport through the AUG H-mode pedestal with applied 3D fields. Our results show that going into the ELM-suppressed regime due to RMPs, impurity transport in the pedestal changes significantly compared to axisymmetric expectations. Impurity diffusivity increases by a factor of 2 in the pedestal. There are radial convection stagnation points both at the bottom and top of the pedestal. Next to enhanced pedestal permeability, W gets trapped in 3D electric potential wells.

Let us first discuss the general validity of our simulations. It has been shown²⁴ that free-boundary JOREK simulations capture well the plasma response and have good agreement with the experiment. The adopted collision operator captures neoclassical impurity transport for the moderate collisionalities of the background plasma. These conditions are satisfied in the pedestal region. No assumptions are made on the W charge distributions as all charge states are modeled and their evolution is governed by effective rate coefficients²⁵ of ionization and recombination. As transport is calculated on a virtually axisymmetric flux surface coordinate, the flux surface corrugation can artificially broaden the density profiles, while particles only follow a deformed flux surface. This however produces only a small effect on the perceived diffusion during simulation startup and does not influence the axisymmetric component of the radial transport.

There are still physics neglected in the model that might alter the transport. As we use a visco-resistive MHD model, it does not include (small scale) turbulence, and effective transport coefficients are applied as closure. Furthermore, a 1-temperature model is deployed for simplicity. Without neutrals, the divertor physics are simplified and JOREK operates in a sheath-limited divertor regime. Further, radiation feedback and electrons seeded by W are not included in this contribution. Finally, W self-collisions are disregarded in the trace impurity limit.

Turbulence will enhance impurity transport. Both due to electric field perturbation and density and temperature fluctuations. Specifically, in the core, it is needed to correctly model transport; however, transport is typically at neoclassical level in the pedestal.⁵² Further, the presence of impurities can cause micro-turbulence stabilization.⁵³ 

For simplicity and to study effects in more isolation, a simplified divertor is used in the simulations. In JOREK, this results in a typically sheath-limited divertor regime. For studying wall-to-core impurity transport in high-recycling regimes, this simplification is likely not sufficient. However, RMP experiments in AUG are performed with low upstream density and rather attached plasma conditions. In this case, the JOREK SOL is similar to the experiment. For obtaining realistic W sources due to sputtering, the divertor geometry and plasma conditions need to be more realistic. JOREK has realistic divertor geometries and kinetic neutrals to obtain realistic divertor conditions.²¹ Next to the source, transport would be altered as the conduction-limited regime is characterized by a parallel temperature gradient in the SOL. This will drive more impurities upstream by the parallel thermal force.

In the trace impurity limit, radiation feedback on the plasma and W self-collisions are neglected. In the used scenarios these assumptions should be acceptable on average. Experiments did not show a decrease in the thermal energy due to tungsten radiation, but locally the plasma could be affected. Local accumulation of W could create pockets of high density where radiation and self-collisions cannot be ignored. However, these two mechanisms will have opposing effects. Radiation cools down the plasma and can deepen the potential well, while self-collisions will limit W accumulation due to enhanced diffusion. Future studies including realistic sputtering sources will need to include at least radiation feedback as impurity densities can locally be much higher in the divertor.

Nonetheless, the current assumptions and limitations, the neglected effects in the studied regions should be small. Further, the main observed transport tendencies remain regardless.

The top of the ITER pedestal will have lower collisionality than AUG, implying our adopted collision operator overestimates the temperature screening at the top of the pedestal. As a result, the flow cancellation of the friction and thermal forces (⁠ $⟨ v neo r ⟩ = 0$⁠) will lie at higher values of $| ∇ T | / T$⁠. To obtain an idea of the implication of using the current collision operator for ITER, let us imagine the normalized gradients—i.e., $| ∇ T | / T$ and $| ∇ n | / n$— of the ITER pedestal are AUG-like as presented in Fig. 8, but with lower collisionality. In Fig. 8, the flow cancelation at the bottom of the pedestal causes radial inward convection to stagnate and accumulate at the bottom of the pedestal. At the top of the pedestal, the flow cancelation point is not an accumulation point. Overestimating H_TSC at the top of the pedestal results in the flow cancelation point lying further inwards. This also means the radial width of effective impurity shielding is larger. If we now assume $H TSC = − 0.3$ in the used AUG pedestal (see Sec. IV) the flow cancelation points change. At the bottom of the pedestal, where stagnation of radial impurity flow causes W accumulation, there is effectively no change. The flow cancelation point at the top of the pedestal is moved radially outwards by about 1 cm in ITER, i.e., $∼ 15$% of the pedestal width. In addition to the neoclassical physics arguments, the 3D ExB vortices arising due to the RMPs should be taken into account in the extrapolation to ITER plasmas, as has been used in the simulations in this work.

Here, we exclude discussing further transport in the core, as more than neoclassical transport needs to be included. More advanced collision operators (such as Ref. 51) include a heat flux limiter that reduces the thermal force and thus the temperature screening effect at low collisionality. This, however, typically includes a free parameter. Further improvement of the collision operator is left for future work.

Once W is in the plasma core it can be measured. The reason our current model cannot compare these results with the experiment rests mostly that turbulent transport is not included. Although the results clearly show transport into the core, we know the transport in that region is not correctly modeled by neoclassical transport alone.

For an indication of pedestal permeability and radial convection speeds through the pedestal, the rise time of a relative radiation signal in the core after LBO could possibly be compared. By comparing when the experimental radiation signal goes up, we might assess average radial convection. It might also be possible to estimate how much of the W reaches the top of the pedestal from the experiment, from this effective shielding could be compared.

Another experimental way to compare is with other impurity species, as lower Z impurity species are observable at lower temperatures. Although we have shown that the Ar-puff does not result in Ar being trapped in potential wells, experimentally measured radial convection through the pedestal can be used to validate the model. Simulations of Ar transport in the pedestal region indicates convection stagnation points lay at the same location, meaning that Ar gathers at the same location in the pedestal. Measuring the location where the Ar transport direction changes sign (the stagnation point) can be used for model validation; and as Ar and W have the stagnation point at the bottom of the pedestal, this can be key for predicting W transport behavior in the ITER pedestal. The difference is that the pedestal is more permeable to Ar. Ar diffuses and transports more quickly into the plasma core. This is consistent with expectations of neoclassical transport. It seems that Ar can thus be experimentally used to study W-like transport through the pedestal with applied 3D fields. Plus, validating our model to experiments will give greater confidence in the future prediction of impurity transport through the pedestal.

This paper introduces the adoption of a collision operator for neoclassical impurity transport valid in the pedestal into the hybrid fluid-kinetic nonlinear MHD code JOREK. We have shown that the collision operator for impurities is valid in the pedestal with moderate values of the background ion collisionality. It correctly simulates neoclassical diffusion and neoclassical convection due to density and temperature gradients for all charge states. The JOREK particle model is now able to simulate collisional neoclassical transport for full-orbit particles, as well as transport due to the electric and magnetic fields. There is no assumption on the existence of flux surfaces, which allows for correct transport simulation on stochastic field lines. Banana–Plateau contributions at low collisionality, Maxwellianization at high collisionality, and i.e., collisions are not included in this model.

With this model, we have been able to investigate tungsten transport in H-mode pedestals with and without applied 3D fields for ELM suppression. The background plasma is modeled after an ASDEX Upgrade (AUG) experiment. By performing a scan over several RMP currents, we have shown the effective impurity diffusion can increase by approximately a factor 2 in the pedestal in the ELM-suppressed regime. This is not the case for RMP currents below the ELM-suppression-threshold, where only a small increase is observed.

Although there is ergodization of the magnetic field in the pedestal, the connection lengths are still $1 – 10$ km. Enhanced poloidal diffusion due to transport along the field line remains limited, due to the high mass of tungsten. Collisional radial convection due to density and temperature gradients is most important in the pedestal. However, the electric field significantly contributes to the radial transport. It shifts the convection stagnation point inward by $∼ 1$ cm onto the separatrix, moving it to a higher temperature.

Depending on the initial W distribution, W is transported more easily through the pedestal, both in- and outwards. For the investigated cases (with relatively deep RMP penetration), RMPs cause a more permeable pedestal for impurity transport as opposed to an axisymmetric pedestal with the same density and temperature profile.

Apart from enhanced diffusion and strong poloidal asymmetries due to RMPs, it causes apparent 3D electric potential wells. Collisional W can get trapped in such potential wells. Notably on the low-field-side outside of the separatrix, close to the X-point. Due to the low temperature, W will not be visible to spectroscopic diagnostic in the region.

To be able to relate closer to the experiments in ASDEX Upgrade a laser-blow-off module has been developed. The quasi-stationary state after the injection results in the same convection-stagnation locations. This results as well in trapping of W outside of the separatrix in an electric potential well. There is no direct experimental validation of W transport through the pedestal with our model. In current experiments, W is only experimentally measured once it is in the plasma core at sufficiently high temperatures.

To investigate whether argon could be used as a substitute for tungsten—as it radiates at lower temperatures—the same type of simulations have been performed for argon with a short puff. Argon does not get trapped in the same potential wells; however, the general transport behavior through the pedestal is similar. Importantly, the radial flow stagnation position at the bottom of the pedestal is the same for argon and tungsten. Furthermore, the pedestal is more permeable to argon and once arrived at the top of the pedestal, radial convection inwards is faster than for tungsten. Our results indicate argon could be used as a substitute for understanding “tungsten-like transport” in the pedestal region with applied 3D fields. Further, dedicated argon experiments could be used to validate our impurity transport model.

Not all facets of impurity transport with applied 3D fields have been discussed. For instance, variation of RMP penetration has not been thoroughly investigated. Future studies will involve a broader scope of different RMP scenarios. Further, simulating and verifying impurity transport in AUG is a stepping stone to predicting tungsten transport in ITER H-mode with applied 3D fields. The current developed model can also be run with radiation feedback on the background plasma. Future studies for ITER can also incorporate heat flux limiters in the collision operator (such as Ref. 51) to account for the reduction of the temperature screening coefficient at lower collisionalities.

To obtain a full impurity transport picture, modeling from wall-to-core is needed, including sources. Realistic divertor conditions are important to predict impurity sputtering and transport, and to study the consequences of transient MHD activity on the divertor. Recent advances³ allow for kinetic neutrals simulation in JOREK to obtain more realistic divertor conditions. In the particle framework, both coupled neutrals and impurities can be included to obtain realistic impurity sources and transport from the divertor to the core.

This work was part of the research programme (No. 19KE01) of the Foundation for Nederlandse Wetenschappelijk Onderzoek Instituten (NWO-I), which is part of the NWO and has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200–EUROfusion) and by the ITER Organization under Implementing Agreement No. 43-2174 with the Eindhoven University of Technology (TU/e) and with TU/e work order No. 107640. ITER is the Nuclear Facility, INB No. 174. Computations were performed using computer clusters Marconi-Fusion and the ITER Scientific Data and Computing Center (SDCC). The views and opinions expressed herein do not necessarily reflect those of the European Commission nor the ITER Organization nor NWO-I. Fruitful discussions with members of IPP Garching (W. Suttrop, R. Dux, and D. Fajardo) are acknowledged. See Ref. 2 for a full list of ASDEX Upgrade team members.

The authors have no conflicts to disclose.

S. Q. Korving: Conceptualization (equal); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (lead). V. Mitterauer: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Software (supporting). G. T. A. Huijsmans: Conceptualization (equal); Funding acquisition (lead); Resources (equal); Software (supporting); Supervision (lead); Writing – review & editing (equal). A. Loarte: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). M. Hoelzl: Data curation (equal); Resources (equal); Software (supporting); Writing – review & editing (equal).

Raw generated data were at the ITER Scientific Data and Computing Center (SDCC). Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

To verify the Homma collision operator, a few benchmark tests are proposed in Ref. 30 to test whether the collision operator can correctly capture the thermal forces on a group of particles. For the exact details of the test cases refer to said publication. The prescribed test cases test the mean acceleration due to the diamagnetic thermal force in scenarios with different temperature gradients, and the mean drift velocity due to the temperature screening effect with different temperature gradients.

For the test simulations I and II Cartesian coordinates are used. The magnetic field $B →$ points to the Z directions. The perpendicular temperature gradient $∇ ⊥$ T is taken along the X-axis. For simplicity the Coulomb logarithm $l n Λ = 15$⁠. Note that the latter is not assumed in the actual simulations.

Test case I is on a short timescale $0 ≤ Ω t ≤ 0.01$⁠. The simulated time evolution of the test particle velocity is directly compared with the theoretical values of the acceleration due to the thermal force. Cases I-1 and I-2 have $∇ ⊥ T = 100$ and 300 eV/m. Case I-3 is the same as I-2 but has an added $∇ ∥ T = 5$ eV/m to simulate simultaneously the diamagnetic and parallel thermal force. Figure 20 shows the average perpendicular velocity for the three different temperature gradients compared to theoretical values. The theoretical values for cases I-2 and I-3 are the same. The JOREK test particles show good agreement with the theory. Figure 21 shows the average parallel velocity over time for case I-3 compared to the theoretical value. The JOREK test particles show good agreement with the theory.

Test case II is on a timescale of many gyrations $0 ≤ Ω t ≤ 700$⁠. The simulated time evolution of the test particles' average position is directly compared with the theoretical values to see whether the model is able to simulate the temperature screening effect.

From Ref. 30, the theoretical velocity of guiding center drift are as follows:

• Case II-1:

  $v scr = − 1.33 e x ( m / s ) ( using Braginskii ) , v scr = − 1.78 e x ( m / s ) ( using Trubnikov ) .$

  (A2)
• Case II-2:

  $v scr = − 4.04 e x ( m / s ) ( using Braginskii ) , v scr = − 5.39 e x ( m / s ) ( using Trubnikov ) .$

  (A3)

Where Braginskii et Trubnikov constitute different derivations of the impurity-background ion collision time τ_ba.

Figure 22 shows the avarage position of the JOREK test particles over time in the X and Y direction for test case II-1. The theoretical Braginskii and Trubnikov values are only for the drift in the X-direction for $∇ ⊥ T = 100$ eV/m. The JOREK test particles show good agreement with theory.

Figure 23 shows the average position of the JOREK test particles over time in the X and Y direction for case II-2. The theoretical Braginskii and Trubnikov values are only for the drift in the X-direction for $∇ ⊥ T = 300$ eV/m. The JOREK test particles show good agreement with theory. The test particles seem to follow more the prediction of the collision time of Braginskii.

The implemented binary collision model (BCM) of Ref. 30 into the JOREK kinetic particle framework has been benchmarked with the prescribed tests. The first test shows the time evolution of the average velocity of the JOREK test particles in the presence of a perpendicular temperature gradient. The results agree with the theoretical values of the thermal force acceleration from kinetic theory. The second test shows the temperature screening effect caused by the thermal force due to a perpendicular temperature gradient. The results agree with theoretical values.

Good agreement of these test simulations shows that JOREK with the added BCM can simulate the thermal force correctly and is reliable enough to be applied for realistic transport simulations.

This appendix supplements the results shown in Sec. III A. Each simulation is done separately, but the results are combined in one plot. The JOREK tungsten particles are initialized with Gaussian profiles with their means at locations $ψ n = 0.28 , 0.5 , 0.72$⁠, and 0.85 and standard deviation $σ = 0.03 Δ ψ n$⁠. The evolution of the density profiles over time as a function of the major radius is shown in Fig. 24. For the first 3 simulations from the left, the means of the distribution are constant. For the most right profile, the means seem to move inwards. This is due to that a part of the particle distribution leaves the computational domain. When this becomes significant, it will distort the estimation of the diffusion coefficient. For Fig. 24, the standard deviations σ can be calculated over time. From this a fit can be made, $σ = 2 D ( t − t 0 )$⁠, to estimate the diffusion coefficient as shown in Fig. 25. The fit in made on partial data. As we initialize the particle profiles with sampling such that T_particle = T_background, the particle distribution is not in an mhd equilibrium with the background plasma. This causes the onset of the expected diffusion to be later, after around $5 − 10$ ms. At the end of the simulation, many particles are lost from the simulation domain. This also distorts the calculation of the profile width. For these reasons, the fit has been done on data between $15 − 45$ ms. Changing the interval slightly does not change the result. The results of Fig. 25 and the comparison with neoclassical theory are presented in Sec. III A.

This appendix section includes supplementary material for the inward pinch test in Sec. III A 2. That subsection presents the average drift of convective velocity of impurities as function of the impurity charge state. These data are derived from the average radial location of the particles of time shown in Fig. 26.