Accumulation of heavy impurities in the tokamak core is detrimental for its performance and can lead to disruption of the plasma. In small to medium-sized tokamaks the effective neoclassical transport in the pedestal is typically oriented radially inward. In larger tokamaks—e.g., ITER—where the temperature gradient is higher and density gradients are lower due to the need to operate in a radiative divertor regime, the neoclassical transport is predicted to be outwards. The models are derived for axisymmetric quasi-steady-state plasmas. Applied 3D magnetic fields, i.e., Resonant Magnetic Perturbation (RMPs) as they are used to suppress Edge Localized Modes (ELMs), have experimentally been observed in AUG to enhance the outflow of heavy impurities in the pedestal. There is no model that can predict neoclassical heavy impurity transport in these ergodized 3D magnetic fields self-consistently. In this contribution, we present our kinetic tungsten transport simulation for an ASDEX Upgrade plasma with applied RMPs. Our model based on Hoelzl et al. [Nucl. Fusion 61, 065001 (2021)], van Vugt et al. [Phys. Plasmas 26, 042508 (2019)], and Korving et al. [Phys. Plasmas 30, 042509 (2023)] utilizes a full-orbit pusher, ionization, recombination, effective line, and continuum radiation and neoclassical collisions with the background plasma. The effective collisional radiative rates are from the OpenADAS database, the neoclassical collision operator uses the framework of Homma et al. [J. Comput. Phys. 250, 206–223 (2013)] and Homma et al. [Nucl. Fusion 56, 036009 (2016)]. We show that the adopted collision operator produces neoclassical transport within a satisfactory degree of accuracy. A sufficiently high RMP current causes an increase in tungsten diffusion in the pedestal by a factor of 2. We compare the average radial transport between axisymmetric and 3D RMP scenarios in the pedestal region. RMPs enhance the pedestal permeability for impurities, which results in enhanced transport. In addition to the enhanced transport, some of W is found to be trapped in 3D potential wells in the scrape-off layer. Due to the lack of suitable diagnostics for W in the pedestal, we investigate and suggest that argon can be an adequate substitute in experiments for model validation and further understanding impurity transport in scenarios with applied 3D magnetic fields. With the newly developed neutral model [Korving et al., Phys. Plasmas 30, 042509 (2023)], we can combine the interaction in the divertor with the 3D RMPs to model the tungsten transport from the divertor toward the core of the plasma.

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 studied1,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.1 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 ELMs1,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 ELMs18 and ELM control by RMPs19 and by pellet triggering.20 JOREK has recently been extended with a kinetic neutral model21 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.

In order to reduce computational requirements, JOREK can also employ a reduced MHD model, which eliminates fast magnetosonic waves while retaining the relevant physics. Reduced MHD has fewer unknowns compared to full MHD, which decreases the computational costs and memory requirements for simulations. The reduced MHD model used in JOREK is derived following the ansatz-based approach.13 The ansatz considered here assumes that the time-dependent part of the magnetic potential is dominated by the toroidal component. The ansatz for the magnetic field is deduced by approximating B ϕ as the vacuum F 0 / R toroidal field. After some calculations we can obtain an expression of the poloidal component of the velocity ( v pol = ( e ϕ × v ) × e ϕ), which captures E × B flows and the ion diamagnetic drift velocity; and the parallel component of the velocity field v = v · B,
B = F 0 R e ϕ + 1 R ψ × e ϕ A · e ϕ = ψ ϕ , v = v pol + v = v E × B + v dia , i + v , = R u × e ϕ ( δ * R / ρ ) p i × e ϕ + v B .
(1)

Here, ψ is the poloidal magnetic flux. The velocity stream function u is defined as Φ / F 0, the electric potential Φ divided by F0. 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.23 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.25 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.

We use a Monte Carlo numerical implementation of the following collision operator in the Landau Form:28,29
( f z t ) coll = b v j e z 2 e b 2 8 π ε 0 2 m z 2 d v [ δ j k u u j u k u 3 ] × [ f z m b f b ( v ) v k f b ( v ) m z f z v k ] .
(2)
With the assumed velocity distribution of the background plasma particles,
f b ( v b ) = n b ( m b 2 π k T b ) 3 / 2 exp ( m b w 2 2 k T b ) × [ 1 m b n b ( k T b ) 2 ( 1 w 2 5 v t h , b 2 ) ( q b heat · w ) ] ,
(3)
where v b is the velocity of a background plasma ion, w v b v b ¯. Note the distribution is a function of v b and not w, thus the Maxwellian is shifted by the local macroscopic velocity v ¯ b. The local background heat flux q b heat distorts the Maxwellian. The two deviations from a Maxwellian result in friction forces (mediated by Coulomb collisions), which lead to neoclassical cross field convective transport.

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.

TABLE I.

Circular plasma parameters used for the collisional neoclassical transport benchmark. The parameters are F0 axial toroidal magnetic field times the major radius, R0 major radius, a minor radius, ne electron density, Te electron temperature, I plasma current. The second table presents the initial tungsten profile position in normalized poloidal flux ψn and the impurity conditionality, where ε 3 / 2 ν z *  1 means the impurity species is in the Pfrisch–Schlüter regime.

F0 (Tm) R0 ( m) a ( m) ne (m−3) Te ( eV) I ( MA)
1020  250  −2.53 
ψ n  0.28  0.5  0.7  0.85   
ε 3 / 2 ν z *   12  14  17  21   
ε 3 / 2 ν i *   0.32  0.38  0.47  0.56   
F0 (Tm) R0 ( m) a ( m) ne (m−3) Te ( eV) I ( MA)
1020  250  −2.53 
ψ n  0.28  0.5  0.7  0.85   
ε 3 / 2 ν z *   12  14  17  21   
ε 3 / 2 ν i *   0.32  0.38  0.47  0.56   

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.31 

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 tgyro and the binary collisions tBCM 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. tcoronal 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 tgyro 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.

FIG. 1.

Advised time steps for resolving gyro-orbits (tgyro) and binary collisions (tBCM) as a function of background plasma temperature over all charge states of W ions with a plasma background density n 0 = 10 20. The two red points indicate the advised time step for tgyro and tBCM at T = 250 eV for the most probable charge state 19. tcoronal is the effective time step needed where W is in coronal equilibrium with the background plasma. The right axis shows the electron and tungsten normalized collisionality.

FIG. 1.

Advised time steps for resolving gyro-orbits (tgyro) and binary collisions (tBCM) as a function of background plasma temperature over all charge states of W ions with a plasma background density n 0 = 10 20. The two red points indicate the advised time step for tgyro and tBCM at T = 250 eV for the most probable charge state 19. tcoronal is the effective time step needed where W is in coronal equilibrium with the background plasma. The right axis shows the electron and tungsten normalized collisionality.

Close modal

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 ne in steady-state.

1. Estimation of neoclassical diffusion

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.

In cylindrical geometry one would expect diffusion to scale as ( Δ X ) 2 / Δ t.5,32 Choosing appropriate step lengths results in the classical diffusion coefficient
D z i C L = ν z i ρ z 2 .
(4)
With ν z i 2 12 π 3 / 2 l n Λ e 4 Z 2 n i ε 0 2 m Z m i v T i 3 the impurity-ion collision frequency,33  ρ z m z v T z e z B the impurity gyro radius and z denoting the impurity species. Including the magnetic topology in a toroidal geometry leads to drift-orbits where particles drift back- and forth around a flux surface. This results in the collisional neoclassical “Pfirsch–Schlüter” diffusion.
D z i P S = q 2 ν z i ρ z 2 ,
(5)
with q, the safety factor.
In the collisional regime, i.e., above the PS-limit ε 3 / 2 ν z * 1 and with only one plasma and impurity species, the neoclassical diffusion becomes the sum of each component3,5,7,8
D z neo = x = C L , B P , P S b a D b a x D z i C L + D z i P S ( 1 + q 2 ) ν z i ρ z 2 ,
(6)
where x is the contribution of the different collisionality regimes—i.e., Classical, Banana-Plateau, Pfirsch–Schlüter— and for each regime the contribution of the other plasma species, a, is summed. The contribution of BP is negligible in the high-collisionality regime. Collisions with electrons are neglected here. Due to the newly introduced collision operator, we expect our diffusion to scale as D z neo / D z C L = 1 + q 2.32,34

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.

FIG. 2.

Diffusion coefficient as a function of the number of collisions per kinetic time step ( t gyro t BCM ) for the benchmark plasma parameters (Table I). The red dotted line indicates the collision amount advised4 to correctly model collisional transport. The shaded areas indicate the relative error to the theoretical diffusion coefficient.

FIG. 2.

Diffusion coefficient as a function of the number of collisions per kinetic time step ( t gyro t BCM ) for the benchmark plasma parameters (Table I). The red dotted line indicates the collision amount advised4 to correctly model collisional transport. The shaded areas indicate the relative error to the theoretical diffusion coefficient.

Close modal

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.

FIG. 3.

Normalized diffusion coefficient as a function of the safety factor. The error bars on the JOREK diffusion coefficient indicate the full-width half-maximum of the initial Gaussian particle distribution. The fit shows JOREK results are marginally higher than the expected neoclassical diffusion coefficient “CL+PS” (Classical + Pfirsch–Schlüter).

FIG. 3.

Normalized diffusion coefficient as a function of the safety factor. The error bars on the JOREK diffusion coefficient indicate the full-width half-maximum of the initial Gaussian particle distribution. The fit shows JOREK results are marginally higher than the expected neoclassical diffusion coefficient “CL+PS” (Classical + Pfirsch–Schlüter).

Close modal

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.

2. Estimation of neoclassical inward pinch and temperature screening effect

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.4 

To simplify and easier test our model with neoclassical theory, the gradient and flows are added artificially on top of the background plasma. In the fluid approximation, collisional neoclassical theory says that in the impurity trace limit the expected radial impurity drift is5,31,35
v neo r 2 q 2 D C L Z z Z i ( K 1 n d n d r H 1 T d T d r ) ,
(7)
where v neo r is the average radial drift, DCL is the neoclassical diffusion coefficient, Zz, Zi is the impurity and background charge state, K is the convective coefficient, H is the impurity screening factor, n, and T is the background plasma density and temperature. In the impurity trace limit, the impurity self-collisions vanish and have no effect on the radial drift. In this subsection, the convective coefficients K = 1, and H = 1 2 are assumed for simplicity of comparison. These values correspond to the simplified large aspect ratio approximation of Pfirsch–Schlüter.36,37
To correctly simulate this radial drift, only including a density and temperature gradient is not enough. For the inward pinch (IWP) we need to include the background-ion fluid velocities arising due to the density gradient. The density gradient drives an ion diamagnetic flow with perpendicular velocity,
v , dia = B × p e n B 2 ,
(8)
where the radial density gradient is included in the pressure gradient p. The diamagnetic flow is not divergence-free. This means another flow must arise to satisfy particle conservation. The so-called Pfirsch–Schlüter (PS) flow in the parallel direction is driven by the diamagnetic flow and satisfies
· ( n i v , P S ) = · ( n i v , dia ) .
(9)
To test the IWP in our model, the background diamagnetic and PS flows are added artificially before sampling the background plasma for the BCM. The velocities for a fixed radial density gradient are then Eq. (8) and
v , P S = F 0 n i e i B p i ψ ( 1 B 2 B 2 ) .
(10)
With F 0 = R 0 B 0, ei, ion background charge, p i ψ, perpendicular ion pressure gradient and B 2 the flux surface average of B2.

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.

FIG. 4.

Average radial velocity due to the inward pinch as a function of the impurity charge state. The black lines indicate the mean of the initial particle distribution, the location of the Half-Maximum of the distribution on the minimum radius (half min) and maximum radius (half max). The solid blue line is the linear fit of the average radial velocity of the particles as a function of the charge state.

FIG. 4.

Average radial velocity due to the inward pinch as a function of the impurity charge state. The black lines indicate the mean of the initial particle distribution, the location of the Half-Maximum of the distribution on the minimum radius (half min) and maximum radius (half max). The solid blue line is the linear fit of the average radial velocity of the particles as a function of the charge state.

Close modal

Figure 4 shows that the JOREK results are satisfactorily close to the neoclassical theory with approximately linear scaling with the impurity charge state Zz. On the low Z side aligns better with the half-max line where qsafe is higher. Zz = 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.

Similarly to the diamagnetic drift, the temperature gradient [from Eq. (7)] drives a diamagnetic heat flux
q , dia = 5 p i 2 e i B 2 B × T ,
(11)
which is not divergence-free in a toroidal geometry and, thus, drives a parallel heat flux,
· ( q , P S ) = · ( q , dia ) .
(12)
This results in the parallel Pfirsch–Schlüter heat flux,
q , P S = 5 F 0 p i 2 e i B T i ψ ( 1 B 2 B 2 ) .
(13)
The thermal force due to this heat flux is what causes the temperature screening effect.

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.

FIG. 5.

Average radial velocity due to the temperature screening effect as a function of the impurity charge state. The black lines indicate the mean of the initial particle distribution, the location of the half-maximum of the distribution on the minimum radius (half min) and maximum radius (half max). The solid blue line is the linear fit of the average radial velocity of the particles as a function of the charge state.

FIG. 5.

Average radial velocity due to the temperature screening effect as a function of the impurity charge state. The black lines indicate the mean of the initial particle distribution, the location of the half-maximum of the distribution on the minimum radius (half min) and maximum radius (half max). The solid blue line is the linear fit of the average radial velocity of the particles as a function of the charge state.

Close modal

3. Effective temperature screening coefficient

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 code38 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-negligible37 and at low collisionalities, there is a reduction of the thermal force due to Banana-Plateau (BP) contributions;37 these result in stronger dependencies on g. FACIT,37 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 HTSC 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 HTSC, FACIT results are shown with the average charge state of the coronal equilibrium at each collisionality.

FIG. 6.

Temperature screening coefficient as a function of the normalized (background) ion collisionality g with constant pressure. The shaded area indicates the collisionality range of the pedestal in AUG during RMP operation. The dashed lines are theoretical models, the fixed large aspect ratio Pfirsch–Schlüter (PS)-model. The Hirshman–Sigmar–Wenzel (HSW)-model includes a fitted dependence on g for impurity-ion collisions. FACIT37 (including Banana-plateau, classical, Pfirsch–Schlüter, and impurity-electron collision terms) with coronal average charge state for W from assumed electron temperature scaling with constant pressure. Approximate most probable coronal charge states are indicated with arrows.

FIG. 6.

Temperature screening coefficient as a function of the normalized (background) ion collisionality g with constant pressure. The shaded area indicates the collisionality range of the pedestal in AUG during RMP operation. The dashed lines are theoretical models, the fixed large aspect ratio Pfirsch–Schlüter (PS)-model. The Hirshman–Sigmar–Wenzel (HSW)-model includes a fitted dependence on g for impurity-ion collisions. FACIT37 (including Banana-plateau, classical, Pfirsch–Schlüter, and impurity-electron collision terms) with coronal average charge state for W from assumed electron temperature scaling with constant pressure. Approximate most probable coronal charge states are indicated with arrows.

Close modal

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,35 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.

TABLE II.

Experimental parameters for ASDEX Upgrade scenario. Time is when the desired ELM suppression occurred. Btor is the nominal toroidal field. Ip is the total plasma current, q95 the safety factor at 95% flux surface, Paux the auxiliary heating (NBI+ECRH), Irmp the RMP current.

Shot No. Time (s) Btor (T) Ip (MA) q95 ( ) Paux ( MW) Irmp ( kA)
#34 548  3.8–3.9  −1.83  0.92  3.61  7.8  1.0 
#39 428  6.3–7.2  −1.83  0.92  3.61  7.8  0.9 
#39 428  4.4–5.1  −1.83  0.92  3.61  7.8  1.01 
#39 428  3.4–4.2  −1.83  0.92  3.61  7.8  1.1 
#39 428  2.5–3.0  −1.83  0.92  3.61  7.8  1.22 
Shot No. Time (s) Btor (T) Ip (MA) q95 ( ) Paux ( MW) Irmp ( kA)
#34 548  3.8–3.9  −1.83  0.92  3.61  7.8  1.0 
#39 428  6.3–7.2  −1.83  0.92  3.61  7.8  0.9 
#39 428  4.4–5.1  −1.83  0.92  3.61  7.8  1.01 
#39 428  3.4–4.2  −1.83  0.92  3.61  7.8  1.1 
#39 428  2.5–3.0  −1.83  0.92  3.61  7.8  1.22 

For ELM suppression in AUG, typically the current in the RMP coils Irmp 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).39 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.24 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.40 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.41 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.42 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

FIG. 7.

Poincare plot on a poloidal cross section of the electron density (a.u.). The toroidal surface indicates variation in the poloidal magnetic flux (a.u.). RMPs create seeds for the islands to form; further, it results in the ergodized layer, and the several plasma lobes below the X-point. There is a kink-response of the plasma, however, not clearly visible in this figure.

FIG. 7.

Poincare plot on a poloidal cross section of the electron density (a.u.). The toroidal surface indicates variation in the poloidal magnetic flux (a.u.). RMPs create seeds for the islands to form; further, it results in the ergodized layer, and the several plasma lobes below the X-point. There is a kink-response of the plasma, however, not clearly visible in this figure.

Close modal

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 IRMP = 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.

FIG. 8.

Profiles characterizing the pedestal on the outer-mid-plane. r l n ( n e ) and | H | r l n ( T e ) indicate the component of perpendicular convective transport. With assumed | H | = 1 / 2 to compare the relative importance of the two components. The arrows indicate the effective convection direction. Locations 1 and 2 indicate where W is initialized in the following simulations. Three regions are indicated for clarity. The Scrape-off-Layer (SOL) outside the separatrix. The pedestal, and the core region. Due to the importance of turbulence in the core, realistic/experimental predictions cannot be made in this region with the current model. The blue/red arrows show the perpendicular convection direction due to collisional neoclassical transport determined by r l n ( n e )- H r l n ( T e ).

FIG. 8.

Profiles characterizing the pedestal on the outer-mid-plane. r l n ( n e ) and | H | r l n ( T e ) indicate the component of perpendicular convective transport. With assumed | H | = 1 / 2 to compare the relative importance of the two components. The arrows indicate the effective convection direction. Locations 1 and 2 indicate where W is initialized in the following simulations. Three regions are indicated for clarity. The Scrape-off-Layer (SOL) outside the separatrix. The pedestal, and the core region. Due to the importance of turbulence in the core, realistic/experimental predictions cannot be made in this region with the current model. The blue/red arrows show the perpendicular convection direction due to collisional neoclassical transport determined by r l n ( n e )- H r l n ( T e ).

Close modal

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, IRMP = 0 and 0.5 kA, where IRMP = 0 is the normal axisymmetric case; and two ELM-suppressed cases with IRMP = 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.

FIG. 9.

(Main) Effective diffusivity of W in the pedestal as a function of the current in the RMP coils. The solid blue line is the expected changed diffusivity due to a change in the ne,Te-profiles while disregarding other changes. The shaded area indicates from which RMP current, ELMs are suppressed. (Insets) Each inset Poincaré plot corresponds to a data point. θpol is the poloidal angle; ψN the normalized flux coordinate with the separatrix at ψ N = 1. The Gaussian particle distributions are initialized at the black dashed line with σ S D = 0.012 [ ψ N ].

FIG. 9.

(Main) Effective diffusivity of W in the pedestal as a function of the current in the RMP coils. The solid blue line is the expected changed diffusivity due to a change in the ne,Te-profiles while disregarding other changes. The shaded area indicates from which RMP current, ELMs are suppressed. (Insets) Each inset Poincaré plot corresponds to a data point. θpol is the poloidal angle; ψN the normalized flux coordinate with the separatrix at ψ N = 1. The Gaussian particle distributions are initialized at the black dashed line with σ S D = 0.012 [ ψ N ].

Close modal

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 IRMP = 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 IRMP = 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.

FIG. 10.

Tungsten transport in the pedestal with IRMP = 1 kA (Top) Time slice of normalized W density at t = 24 ms. (Bottom) flux-surface-averaged normalized W density in the pedestal evolving over time. (a) simulation including the effects of the E-field and collisions on the transport. (b) Excluding E-field. (c) Excluding the collisions. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

FIG. 10.

Tungsten transport in the pedestal with IRMP = 1 kA (Top) Time slice of normalized W density at t = 24 ms. (Bottom) flux-surface-averaged normalized W density in the pedestal evolving over time. (a) simulation including the effects of the E-field and collisions on the transport. (b) Excluding E-field. (c) Excluding the collisions. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

Close modal

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, vr, 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.

FIG. 11.

Contours of electric potential and direction of E × B drift in the edge region on top of the electron density. White contours indicate the potential contours without the application of RMPS. Black contours indicate the perturbed potential profiles. The pink contour is the unperturbed separatrix. The arising E × B vortices all rotate in the same direction.

FIG. 11.

Contours of electric potential and direction of E × B drift in the edge region on top of the electron density. White contours indicate the potential contours without the application of RMPS. Black contours indicate the perturbed potential profiles. The pink contour is the unperturbed separatrix. The arising E × B vortices all rotate in the same direction.

Close modal

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.

FIG. 12.

Tungsten transport at the bottom of the pedestal with (left) and without (right) RMPs at IRMP = 1 kA (top) time slice of normalized W density at t = 22.4 ms. (Bottom) flux-surface-averaged normalized W density evolving over time. Initialized just outside of the separatrix (a) simulation included the effects of RMPs (b) Axisymmetric case without RMPs. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

FIG. 12.

Tungsten transport at the bottom of the pedestal with (left) and without (right) RMPs at IRMP = 1 kA (top) time slice of normalized W density at t = 22.4 ms. (Bottom) flux-surface-averaged normalized W density evolving over time. Initialized just outside of the separatrix (a) simulation included the effects of RMPs (b) Axisymmetric case without RMPs. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

Close modal

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.

FIG. 13.

Tungsten transport in the pedestal with (left) and without (right) RMPs at IRMP = 1 kA (Top) Time slice of normalized W density at t = 8.5 ms. (Bottom) Flux-surface-averaged normalized W density evolving over time. Initialized just inside of the separatrix (a) simulation included the effects of RMPs (b) Axisymmetric case without RMPs. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

FIG. 13.

Tungsten transport in the pedestal with (left) and without (right) RMPs at IRMP = 1 kA (Top) Time slice of normalized W density at t = 8.5 ms. (Bottom) Flux-surface-averaged normalized W density evolving over time. Initialized just inside of the separatrix (a) simulation included the effects of RMPs (b) Axisymmetric case without RMPs. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

Close modal

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.31 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 IRMP = 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.

FIG. 14.

Poloidal tungsten density cross sections and contours of tungsten density during a quasi-steady-state part of a simulation. Solid Pink: highest W density. Transparent light pink: relatively high W density. Due to the applied 3D fields, W mostly tends to gather at the top and high-field side. However, a pronounced toroidally localized trapping takes place at the bottom of the pedestal right of the X-point.

FIG. 14.

Poloidal tungsten density cross sections and contours of tungsten density during a quasi-steady-state part of a simulation. Solid Pink: highest W density. Transparent light pink: relatively high W density. Due to the applied 3D fields, W mostly tends to gather at the top and high-field side. However, a pronounced toroidally localized trapping takes place at the bottom of the pedestal right of the X-point.

Close modal

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.

FIG. 15.

Electric potential (left axis) and projection of W density (right axis) along the toroidal direction ϕ, for a location where it intersects with the trapped W blob (R,Z = 1.567 m, −0.788 m, about 20 cm top-right from the unperturbed X-point). Tungsten gathers in the potential well. Vertical axes in au with different zero-lines. Oscillations around zero in the W density are an artifact due to the projection not fully resolving the sharp gradients.

FIG. 15.

Electric potential (left axis) and projection of W density (right axis) along the toroidal direction ϕ, for a location where it intersects with the trapped W blob (R,Z = 1.567 m, −0.788 m, about 20 cm top-right from the unperturbed X-point). Tungsten gathers in the potential well. Vertical axes in au with different zero-lines. Oscillations around zero in the W density are an artifact due to the projection not fully resolving the sharp gradients.

Close modal

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.47 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 IRMP = 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).48 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.

FIG. 16.

n = 0 component of W density projections of the initial phase of a W-LBO simulation. The “cold” W speeds up to a convection quasi-stationary state in about 7 ms. On the longer scale diffusion determines the transport. The color bar is logarithmic. The dottedness of W in the core region is a result of a finite particle amount and logarithmic coloring. Oscillations and line artifacts are due to the projection method and do not reflect the actual location of W.

FIG. 16.

n = 0 component of W density projections of the initial phase of a W-LBO simulation. The “cold” W speeds up to a convection quasi-stationary state in about 7 ms. On the longer scale diffusion determines the transport. The color bar is logarithmic. The dottedness of W in the core region is a result of a finite particle amount and logarithmic coloring. Oscillations and line artifacts are due to the projection method and do not reflect the actual location of W.

Close modal

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.

FIG. 17.

3D rendering of poloidal cross sections of the background electron density (front) and of the W density in the SOL and edge (back). The gray contours show how the bulk of the W is spread at around t t 0 = 0.5 ms.

FIG. 17.

3D rendering of poloidal cross sections of the background electron density (front) and of the W density in the SOL and edge (back). The gray contours show how the bulk of the W is spread at around t t 0 = 0.5 ms.

Close modal

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., IRMP, 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.

FIG. 18.

n = 0 projection of Ar density sufficiently long after Ar-puff to obtain a quasi-stationary state. The dottedness of W in the core region is a result of a finite particle number and logarithmic coloring. Oscillations and line artifacts are due to the projection method and do not reflect the actual location of Ar.

FIG. 18.

n = 0 projection of Ar density sufficiently long after Ar-puff to obtain a quasi-stationary state. The dottedness of W in the core region is a result of a finite particle number and logarithmic coloring. Oscillations and line artifacts are due to the projection method and do not reflect the actual location of Ar.

Close modal

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.

FIG. 19.

Argon transport in the pedestal with RMPs at IRMP = 1 kA (Top) Time slice of normalized Ar density at t = 28.2 ms. (Bottom) Flux-surface-averaged normalized Ar density evolving over time. (a) Initialized just inside of the separatrix. (b) Initialized just outisde of the separatrix. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

FIG. 19.

Argon transport in the pedestal with RMPs at IRMP = 1 kA (Top) Time slice of normalized Ar density at t = 28.2 ms. (Bottom) Flux-surface-averaged normalized Ar density evolving over time. (a) Initialized just inside of the separatrix. (b) Initialized just outisde of the separatrix. Δ r separatrix is the distance from the (unperturbed) separatrix at the outer mid-plane (OMP).

Close modal

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 shown24 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 coefficients25 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.52 Further, the presence of impurities can cause micro-turbulence stabilization.53 

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.21 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 HTSC 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 advances3 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.

1. Test simulation I—Diamagnetic thermal force

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.

FIG. 20.

Results of test cases (I-1) (I-2), and (I-3) and theoretical values “theo”: The Y-component of the average velocity of JOREK test particles accelerated only by the diamagnetic thermal force v a , Y ¯ T The Y-component of the average velocity of case-3 should and has the same result as case (I-2).

FIG. 20.

Results of test cases (I-1) (I-2), and (I-3) and theoretical values “theo”: The Y-component of the average velocity of JOREK test particles accelerated only by the diamagnetic thermal force v a , Y ¯ T The Y-component of the average velocity of case-3 should and has the same result as case (I-2).

Close modal
FIG. 21.

Result of test case (I-3) and theoretical value “Theo”: The Z-component of the average velocity of JOREK test particles accelerated by the parallel thermal force v a , Z ¯ T.

FIG. 21.

Result of test case (I-3) and theoretical value “Theo”: The Z-component of the average velocity of JOREK test particles accelerated by the parallel thermal force v a , Z ¯ T.

Close modal
2. Test simulation II—Temperature screening effect

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.

Two perpendicular temperature gradients have been used,
( II 1 ) T = 100 eV / m e X , ( II 2 ) T = 300 eV / m e X .
(A1)

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.

FIG. 22.

Result of case (II-1): JOREK X and JOREK Y represent the average value of the simulated X-position X ¯ ( t ) and the Y-position Y ¯ ( t ) of the test particles. The Braginskii estimation using τ b a B and the Trubnikov estimation by τ b a T are plotted, respectively by the dashed and the solid line.

FIG. 22.

Result of case (II-1): JOREK X and JOREK Y represent the average value of the simulated X-position X ¯ ( t ) and the Y-position Y ¯ ( t ) of the test particles. The Braginskii estimation using τ b a B and the Trubnikov estimation by τ b a T are plotted, respectively by the dashed and the solid line.

Close modal

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.

FIG. 23.

Result of case (II-2): Closed squares and cross marks represent the average value of the simulated X-position X ¯ ( t ) and the Y-position Y ¯ ( t ) of the JOREK test particles. The Braginskii estimation using τ b a B and the Trubnikov estimation by τ b a T are plotted, respectively by the dashed and the solid line.

FIG. 23.

Result of case (II-2): Closed squares and cross marks represent the average value of the simulated X-position X ¯ ( t ) and the Y-position Y ¯ ( t ) of the JOREK test particles. The Braginskii estimation using τ b a B and the Trubnikov estimation by τ b a T are plotted, respectively by the dashed and the solid line.

Close modal
3. Summary

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 Tparticle = Tbackground, 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.

FIG. 24.

Evolution of projected particle density as a function of major radius over time for four different initial positions.

FIG. 24.

Evolution of projected particle density as a function of major radius over time for four different initial positions.

Close modal
FIG. 25.

Evolution of standard deviation of the Gaussian density profiles over time. For each simulation, the diffusion coefficient is calculated from a fit of the standard deviations over time between 15 and 45 ms. The exact interval does not change the result, only if the first 10 ms and the last 5 ms are excluded.

FIG. 25.

Evolution of standard deviation of the Gaussian density profiles over time. For each simulation, the diffusion coefficient is calculated from a fit of the standard deviations over time between 15 and 45 ms. The exact interval does not change the result, only if the first 10 ms and the last 5 ms are excluded.

Close modal

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.

FIG. 26.

Average radial position of impurities over time for several different charge states. Due to the background density gradient and resulting diamagnetic and Pfirsch–Schlüter flows the impurities exhibit an inward pinch.

FIG. 26.

Average radial position of impurities over time for several different charge states. Due to the background density gradient and resulting diamagnetic and Pfirsch–Schlüter flows the impurities exhibit an inward pinch.

Close modal
1.
D. C.
van Vugt
,
G. T. A.
Huijsmans
,
M.
Hoelzl
,
A.
Loarte
,
ASDEX Upgrade Team
, and
Eurofusion MST1 Team
, “
Kinetic modeling of ELM-induced tungsten transport in a tokamak plasma
,”
Phys. Plasmas
26
,
042508
(
2019
).
2.
H.
Zohm
et al, “
Overview of ASDEX Upgrade results in view of ITER and DEMO
,”
Nucl. Fusion
(published online) (
2024
).
3.
F. L.
Hinton
and
R. D.
Hazeltine
, “
Theory of plasma transport in toroidal confinement systems
,”
Rev. Mod. Phys.
48
,
239
(
1976
).
4.
S. P.
Hirshman
and
D. J.
Sigmar
, “
Neoclassical transport of impurities in tokamak plasmas
,”
Nucl. Fusion
21
,
1079
(
1981
).
5.
P.
Helander
and
D. J.
Sigmar
,
Collisional Transport in Magnetized Plasmas
(
Cambridge University Press
,
2002
).
6.
V. M.
Zhdanov
,
Transport Processes in Multicomponent Plasma
(
Elsevier Science Publishers B.V
.,
2002
).
7.
R.
Balescu
,
Transport Processes in Plasmas: 2. Neoclassical Transport Theory
(
Elsevier Science Publishers B.V
.,
1988
) Chap. 12, 13, and 15.
8.
R.
Dux
,
A.
Loarte
,
E.
Fable
, and
A.
Kukushkin
, “
Transport of tungsten in the H-mode edge transport barrier of ITER
,”
Plasma Phys. Controlled Fusion
56
,
124003
(
2014
).
9.
F.
Wang
,
X. J.
Zha
,
Y. M.
Duan
,
S. T.
Mao
,
L.
Wang
,
F.
Zhong
,
Y.
Liang
,
L.
Li
,
H. W.
Lu
,
L. Q.
Hu
,
Y. P.
Chen
, and
Z. D.
Yang
, “
Simulations on W impurity transport in the edge of EAST H-mode plasmas
,”
Plasma Phys. Controlled Fusion
60
,
125005
(
2018
).
10.
G.
Vogel
,
H.
Zhang
,
Y.
Shen
,
S.
Dai
,
Y.
Sun
,
J.
Huang
,
S.
Gu
,
J.
Fu
,
R.
Hu
,
J.
Chen
,
X.
Du
,
Q.
Wang
,
Y.
Yu
,
S.
Mao
,
B.
Lyu
, and
M.
Ye
, “
Experimental and simulation study of impurity transport response to RMPs in RF-heated H-mode plasmas at EAST
,”
J. Plasma Phys.
87
,
905870213
(
2021
).
11.
S. Y.
Dai
,
H. M.
Zhang
,
B.
Lyu
,
Y. W.
Sun
,
M. N.
Jia
,
Y.
Feng
,
Z. X.
Wang
, and
D. Z.
Wang
, “
EMC3-EIRENE modelling of tungsten behavior under resonant magnetic perturbations on EAST: Effects of tungsten sputtering and impurity screening
,”
Plasma Phys. Controlled Fusion
63
,
025003
(
2020
).
12.
P.
Sinha
,
N. M.
Ferraro
, and
E.
Belli
, “
Neoclassical transport due to resonant magnetic perturbations in DIII-D
,”
Nucl. Fusion
62
,
126028
(
2022
).
13.
M.
Hoelzl
,
G. T. A.
Huijsmans
,
S. J. P.
Pamela
,
M.
Bécoulet
,
E.
Nardon
,
F. J.
Artola
,
B.
Nkonga
,
C. V.
Atanasiu
,
V.
Bandaru
,
A.
Bhole
et al, “
The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasma
,”
Nucl. Fusion
61
,
065001
(
2021
).
14.
G. T. A.
Huijsmans
and
O.
Czarny
, “
MHD stability in X-point geometry: Simulation of ELMs
,”
Nucl. Fusion
47
,
659
(
2007
).
15.
O.
Czarny
and
G. T. A.
Huijsmans
, “
Bézier surfaces and finite elements for MHD simulations
,”
J. Comput. Phys.
227
,
7423
(
2008
).
16.
M.
Hoelzl
,
A.
Cathey
,
S.
Futatani
,
G. T. A.
Huijsmans
,
F.
Orain
,
M.
Dunne
,
S.
Pamela
,
M.
Becoulet
,
F. J.
Artola
,
D. C.
van Vugt
,
S.
Smith
,
N.
Schwarz
,
F.
Liu
,
S. Q.
Korving
et al, “
Simulations of Edge Localized Mode (ELM) Cycles and ELM Control
,” in
28th IAEA Fusion Energy Conference (FEC)
,
2020
.
17.
G. T. A.
Huijsmans
,
S.
Pamela
,
E.
van der Plas
, and
P.
Ramet
, “
Non-linear MHD simulations of edge localized modes (ELMs)
,”
Plasma Phys. Controlled Fusion
51
,
124012
(
2009
).
18.
A.
Cathey
,
M.
Hoelzl
,
K.
Lackner
,
G. T. A.
Huijsmans
,
M. G.
Dunne
,
E.
Wolfrum
,
S. J. P.
Pamela
,
F.
Orain
,
S.
Günter
,
JOREK Team
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Non-linear extended MHD simulations of type-I edge localised mode cycles in ASDEX Upgrade and their underlying triggering mechanism
,”
Nucl. Fusion
60
,
124007
(
2020
).
19.
M.
Becoulet
,
G. T. A.
Huijsmans
,
C.
Passeron
,
Y. Q.
Liu
,
T. E.
Evans
,
L. L.
Lao
,
L.
Li
,
A.
Loarte
,
S. D.
Pinches
,
A.
Polevoi
,
M.
Hosokawa
,
S. K.
Kim
,
S. J. P.
Pamela
,
S.
Futatani
, and
JOREK Team
, “
Non-linear MHD modelling of edge localized modes suppression by resonant magnetic perturbations in ITER
,”
Nucl. Fusion
62
,
066022
(
2022
).
20.
S.
Futatani
,
A.
Cathey
,
M.
Hoelzl
,
P. T.
Lang
,
G. T. A.
Huijsmans
,
M.
Dunne
,
JOREK Team
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Transition from no-ELM response to pellet ELM triggering during pedestal build-up–insights from extended MHD simulations
,”
Nucl. Fusion
61
,
046043
(
2021
).
21.
S. Q.
Korving
,
G. T. A.
Huijsmans
,
J.-S.
Park
,
A.
Loarte
, and
JOREK Team
, “
Development of the neutral model in the nonlinear MHD code JOREK: Application to E × B drifts in ITER PFPO-1 plasmas featured
,”
Phys. Plasmas
30
,
042509
(
2023
).
22.
D. C.
van Vugt
, Ph.D. thesis, Eindhoven University of Technology,
2019
.
23.
M.
Hölzl
,
P.
Merkel
,
G. T. A.
Huysmans
,
E.
Nardon
,
E.
Strumberger
,
R.
McAdams
,
I.
Chapman
,
S.
Günter
, and
K.
Lackner
, “
Coupling JOREK and STARWALL codes for non-linear resistive-wall simulations
,”
J. Phys.
401
,
012010
(
2012
).
24.
V.
Mitterauer
,
M.
Hoelzl
,
M.
Willensdorfer
,
M.
Dunne
,
N.
Schwarz
,
J.
Artola
,
JOREK Team
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Non-linear free boundary simulations of the plasma response to resonant magnetic perturbations in ASDEX Upgrade plasmas
,”
J. Phys.
2397
,
012008
(
2022
).
25.
H.
Summers
, see http://www.adas.ac.uk for “
The ADAS User Manual
” (
2004
).
26.
H.
Qin
,
S.
Zhang
,
J.
Xiao
,
J.
Liu
,
Y.
Sun
, and
W. M.
Tang
, “
Why is Boris algorithm so good
,”
Phys. Plasmas
20
,
084503
(
2013
).
27.
G. L.
Delzanno
and
E.
Camporeale
, “
On particle movers in cylindrical geometry for Particle-In-Cell simulations
,”
J. Comput. Phys.
253
,
259
277
(
2013
).
28.
T.
Takizuka
and
H.
Abe
, “
A binary collision model for plasma simulation with a particle code
,”
J. Comput. Phys.
25
,
205
219
(
1977
).
29.
Y.
Homma
and
A.
Hatayama
, “
Numerical modeling of thermal force in a plasma for test-ion transport simulation based on Monte Carlo Binary Collision Model
,”
J. Comput. Phys.
231
,
3211
3227
(
2012
).
30.
Y.
Homma
and
A.
Hatayama
, “
Numerical modeling of the thermal force in a plasma for test-ion transport simulation based on a Monte Carlo Binary Collision Model (II)—Thermal forces due to temperature gradients parallel and perpendicular to the magnetic field
,”
J. Comput. Phys.
250
,
206
223
(
2013
).
31.
P.
Maget
,
J.
Frank
,
T.
Nicolas
,
O.
Agullo
,
X.
Garbet
, and
H.
Lütjens
, “
Natural poloidal asymmetry and neoclassical transport of impurities in tokamak plasmas
,”
Plasma Phys. Controlled Fusion
62
,
025001
(
2020
).
32.
K. T.
Tsang
,
Y.
Matsuda
, and
H.
Okuda
, “
Numerical simulation of neoclassical diffusion
,”
Phys. Fluids
18
,
1282
1286
(
1975
).
33.
W.
Fundamenski
and
O. E.
Garcia
, “
Comparison of Coulomb collision rates in the plasma physics and magnetically confined fusion literature
,”
Report No. EFDA-JET-R(07)01
(EFDA-JET,
2007
), https://scipub.euro-fusion.org/archives/jet-archive/comparison-of-coulomb-collision-rates-in-the-plasma-physics-and-magnetically-confined-fusion-literature.
34.
Y.
Sawada
,
M.
Toma
,
Y.
Homma
,
W.
Sato
,
T.
Furuta
,
S.
Yamoto
, and
A.
Hatayama
, “
Modeling of impurity classical/neoclassical transport by Monte-Carlo binary collision model
,”
Fusion Sci. Technol.
63
,
352
354
(
2013
).
35.
Y.
Homma
,
S.
Yamoto
,
Y.
Sawada
,
H.
Inoue
, and
A.
Hatayama
, “
Kinetic modelling for neoclassical transport of high-Z impurity particles using a binary collision method
,”
Nucl. Fusion
56
,
036009
(
2016
).
36.
A.
Samain
and
F.
Werkoff
, “
Diffusion in tokamaks with impurities in the Pfirsch-Schlüter regime
,”
Nucl. Fusion
17
,
53
(
1977
).
37.
D.
Fajardo
,
C.
Angioni
,
P.
Maget
, and
P.
Manas
, “
Analytical model for collisional impurity transport in tokamaks at arbitrary collisionality
,”
Plasma Phys. Controlled Fusion
64
,
055017
(
2022
).
38.
E. A.
Belli
and
J.
Candy
, “
Full linearized Fokker–Planck collisions in neoclassical transport simulations
,”
Plasma Phys. Controlled Fusion
54
,
015015
(
2012
).
39.
W.
Suttrop
,
A.
Kirk
,
V.
Bobkov
,
M.
Cavedon
,
M.
Dunne
,
R. M.
McDermott
,
H.
Meyer
,
R.
Nazikian
,
C.
Paz-Soldan
, and
D. A.
Ryan
, “
Experimental conditions to suppress edge localised modes by magnetic perturbations in the ASDEX Upgrade tokamak
,”
Nucl. Fusion
58
,
096031
(
2018
).
40.
Q. M.
Hu
,
R.
Nazikian
,
B. A.
Grierson
,
N. C.
Logan
,
J.-K.
Park
,
C.
Paz-Soldan
, and
Q.
Yu
, “
The density dependence of edge-localized-mode suppression and pump-out by resonant magnetic perturbations in the DIII-D tokamak
,”
Phys. Plasmas
26
,
120702
(
2019
).
41.
N.
Leuthold
,
W.
Suttrop
,
M.
Willensdorfer
,
G.
Birkenmeier
,
D.
Brida
,
M.
Cavedon
,
M.
Dunne
,
G.
Conway
,
R.
Fischer
,
L.
Gil
,
T.
Happel
,
P.
Hennequin
,
A.
Kappatou
,
A.
Kirk
,
P.
Manz
,
R.
McDermott
,
J.
Vicente
,
H.
Zohm
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Turbulence characterization during the suppression of edge-localized modes by magnetic perturbations on ASDEX upgrade
,”
Nucl. Fusion
63
,
046014
(
2023
).
42.
M.
Pueschel
,
D.
Hatch
,
M.
Kotschenreuther
,
A.
Ishizawa
, and
G.
Merlo
, “
Multi-scale interactions of microtearing turbulence in the tokamak pedestal
,”
Nucl. Fusion
60
,
124005
(
2020
).
43.
M.
Bécoulet
,
F.
Orain
,
G. T. A.
Huijsmans
,
S.
Pamela
,
P.
Cahyna
,
M.
Hoelzl
,
X.
Garbet
,
E.
Franck
,
E.
Sonnendrücker
et al, “
Mechanism of edge localized mode mitigation by resonant magnetic perturbations
,”
Phys. Rev. Lett.
113
,
115001
(
2014
).
44.
F.
Orain
,
M.
Hoelzl
,
F.
Mink
,
M.
Willensdorfer
,
M.
Bécoulet
,
M.
Dunne
,
S.
Günter
,
G.
Huijsmans
,
K.
Lackner
,
S.
Pamela
,
W.
Suttrop
,
E.
Viezzer
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Non-linear modeling of the threshold between ELM mitigation and ELM suppression by resonant magnetic perturbations in ASDEX upgrade
,”
Phys. Plasmas
26
(
4
),
042503
(
2019
).
45.
R.
Fitzpatrick
, “
Bifurcated states of a rotating tokamak plasma in the presence of a static error-field
,”
Phys. Plasmas
5
,
3325
3341
(
1998
).
46.
M.
Willendorfer
,
V.
Mitterauer
,
M.
Hoelzl
,
W.
Suttrop
,
M.
Cianciosa
,
M.
Dunne
,
R.
Fischer
,
N.
Leuthold
,
O.
Samoylov
,
G.
Suarez Lopez
, and
D.
Wendler
, “
Magnetic islands in 3D tokamak plasmas during the suppression of filamentary eruptions
,”
Research Square
preprint (2023).
47.
T.
Lunt
,
Y.
Feng
,
K.
Krieger
,
A.
Kallenbach
,
R.
Neu
,
A.
Janzer
,
R.
Dux
,
D. P.
Coster
,
M.
Wischmeier
,
N.
Hicks
,
W.
Suttrop
,
T.
Pütterich
,
H.
Müller
,
V.
Rohde
,
E.
Wolfrum
,
R.
Fischer
, and
ASDEX Upgrade Team
, “
3D modeling of the ASDEX Upgrade edge plasma exposed to a localized tungsten source by means of EMC3-Eirene
,”
J. Nucl. Mater.
415
,
S505
S508
(
2011
).
48.
F.
Ryter
,
R.
Neu
,
R.
Dux
,
H.-U.
Fahrbach
,
F.
Leuterer
,
G.
Pereverzev
,
J.
Schweinzer
,
J.
Stober
,
W.
Suttrop
,
ASDEX Upgrade Team
,
F. D.
Luca
,
A.
Jacchia
, and
J. E.
Kinsey
, “
Propagation of cold pulses and heat pulses in ASDEX upgrade
,”
Nucl. Fusion
40
,
1917
(
2000
).
49.
D.
Reiser
,
D.
Reiter
, and
M. Z.
Tokar
, “
Improved kinetic test particle model for impurity transport in tokamaks
,”
Nucl. Fusion
38
,
165
(
1998
).
50.
R.
Dux
,
Impurity Transport in Tokamak Plasmas
(
Max-Planck-Institut fur Plasmaphysik
,
2004
).
51.
Y.
Homma
,
K.
Hoshino
,
S.
Yamoto
,
S.
Tokunaga
,
N.
Asakura
,
Y.
Sakamoto
, and
Joint Special Design Team for Fusion DEMO
, “
Effect of collisionality dependence of thermal force on impurity transport under a lower collisional condition in demo scrape-off layer plasma
,”
Nucl. Fusion
60
,
046031
(
2020
).
52.
E.
Viezzer
,
M.
Cavedon
,
E.
Fable
,
F. M.
Laggner
,
R. M.
McDermott
,
J.
Galdon-Quiroga
,
M. G.
Dunne
,
A.
Kappatou
,
C.
Angioni
,
P.
Cano-Megias
,
D. J.
Cruz-Zabala
,
R.
Dux
,
T.
Pütterich
,
F.
Ryter
,
E.
Wolfrum
,
ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Ion heat transport dynamics during edge localized mode cycles at ASDEX Upgrade
,”
Nucl. Fusion
58
,
026031
(
2018
).
53.
M.
Marin
,
J.
Citrin
,
C.
Giroud
,
C.
Bourdelle
,
Y.
Camenen
,
L.
Garzotti
,
A.
Ho
,
M.
Sertoli
, and
JET Contributors
, “
Integrated modelling of neon impact on JET H-mode core plasmas
,”
Nucl. Fusion
63
,
016019
(
2022
).