Effect of toroidal rotation on impurity transport in tokamak improved confinement

The centrifugal force effects from toroidal rotation in improved confinement plasmas are analyzed on high-Z impurities in tokamaks. Tungsten (W) transport simulations are performed using the impurity transport code developed in the integrated code TASK. The geometric factors PA and PB are introduced into the neoclassical transport coefficients to include the effects of the toroidal rotation, which come from poloidal asymmetry in the high-Z impurity profiles. Inward neoclassical particle pinch driven by the main ion density gradient is enhanced by the poloidal asymmetry to be the dominant mechanism for W accumulation in the plasma central region. Simulations with experimental plasma profiles show good agreement with the experimental results and first-principle simulation results in the H-mode. In the hybrid mode and advanced mode, the impurity accumulation is enhanced in the internal transport barrier (ITB) regions. The condition to suppress impurity accumulation is investigated by calculating dependencies on the toroidal rotation velocity and ITB position. The neoclassical transport is sufficiently small with the prospected ITER condition of the Mach number of main ions Mi ∼ 0.1. The impurity transport inside the ITB is strongly influenced by competition between the density peaking effect and the temperature screening effect, and the present simulations show suppression of the impurity accumulation with the outer ITB position to improve the plasma performance, due to the relatively larger temperature gradient of the main ion.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Various improved confinement modes observed in tokamak plasma experiments contribute to steady-state confinement of high-temperature and high-density plasmas.The high confinement mode (H-mode), in which turbulence is suppressed in the edge region and density and temperature are increased in the confinement region, has been adopted as an operating scenario in ITER.Higher plasma performance can be achieved with formation of an internal transport barrier (ITB) in the plasma confinement region, and is being studied as an advanced scenario to achieve steady-state operation in ITER [1][2][3].The hybrid scenario, which is an intermediate operation mode, has been investigated in several experimental devices [4][5][6].On the other hand, significant degradation of core plasma performance has been reported due to radiation losses as a result of accumulation of high-Z impurities such as tungsten (W), which is used for the divertor and vacuum vessel wall material in fusion reactors [7][8][9].In particular, the impurity accumulation can be accelerated in the ITB region of improved confinement discharges in negative magnetic shear mode and high-β p mode.One of the key mechanisms for the impurity accumulation in plasmas is inward neoclassical transport driven by the bulk ion density gradient.Moreover, it has been experimentally and theoretically confirmed that the neoclassical transport is enhanced by the centrifugal effects due to the toroidal rotation of the bulk plasma, which causes poloidal asymmetry of the W distribution [10][11][12][13].The transport of high-Z impurities is dominated by neoclassical transport rather than turbulent transport in the central region of the plasma, and is enhanced by the inward neoclassical pinch in ITB plasmas with suppressed turbulence [14][15][16][17][18][19].An effective method to control the accumulation of impurities is being studied, in which radiofrequency (RF) heating is applied into the core plasma to reduce the neoclassical inward pinch [20][21][22].An integrated analytical method is effective for comprehensive prediction of physical phenomena observed in fusion plasmas.Various integrated codes combining physical and engineering models have been developed, and are useful not only for analyzing experiment results, but also for understanding physical phenomena [23][24][25][26][27]. Plasma collapse due to impurities must be avoided.Understanding and predicting impurity behavior in tokamaks is important to solve the problem of impurity transport.In order to establish operation scenarios including impurity control, a fast and accurate prediction tool for entire plasma discharges is required both for core plasma and impurity behavior in integrated manner.
The purpose of this study is to investigate impurity transport behavior in improved confinement tokamak plasmas forming ITBs.Variation of the neoclassical transport is paid attention, and plays an important role in the plasma central regions.We have developed a 1.5D tokamak transport simulation code including impurity transport in the integrated transport simulation code TASK (Transport Analyzing System for tokamaK) [28,29].Simulations focused on the W neoclassical transport are performed for a tokamak experimental device, corresponding to the three confinement modes considered in the ITER operation scenario.As the neoclassical transport model, the NCLASS code is used, and since it does not include the effects of impurity poloidal asymmetry, the model is extended to include the effects of the toroidal rotation.Then, impurity transport simulations are performed to investigate the effects of the toroidal rotation on impurity distribution formation in improved confinement modes.Magnitudes of impurities accumulation are evaluated for the temperature screening effects in the improvement transport region.

Plasma transport model
The TASK code is one of the integrated codes developed in Japan, and has a modular structure combining many codes such as equilibrium, transport, and wave analyses.The multiple modules exchange data with each other through the data interface BPSD [30], which enables self-consistent simulations of time evolutions of plasma quantities.
In the transport calculations of plasma main component using the TASK/TR code (transport module for main plasma species), the following one-dimensional diffusive equations for particle and heat transport and magnetic diffusion are solved.
∂ ∂t Here, n s , T s , S s , P s , D s , V s , χ s , V Es , B t , η // , J CD , J BS , I, ρ = r/a, are the plasma density, the plasma temperature, the particle source, the heat source, the particle diffusion coefficient, the particle convective velocity, the heat diffusion coefficient, the heat convective velocity, the toroidal magnetic field, the plasma resistivity along the magnetic field lines, the induced current density, the bootstrap current, B t R and the normalized minor radius, respectively, and V' = dV/dρ is the radial derivative of plasma volume and < > denotes a magnetic flux surface average.In the particle transport equation, the quasi-neutrality condition is assumed.Transport coefficients include neoclassical and turbulent contributions.
The heating source of neutral beam injection (NBI) is given by Here, P NB,0 is the intensity of NBI power, ρ NB 0 is the radial position of NBI, ρ NB w is the radial width of NBI.The lower hybrid (LH) heating profile is also given with the same form of the equation.A H-mode pedestal profile is given by a modified tanh function expressed as [31] X Here, X represents each plasma quantity, and X ped , X sol , ρ mid , ρ w and α ped are the height of pedestal, the height of scrapeoff layer, the middle position of pedestal, the radial width of pedestal, the radial profile fitting parameter, respectively.The diffusion coefficients are given by linear summation of neoclassical and turbulent diffusion coefficients.For the neoclassical transport model, NCLASS [32] is applied.For the turbulent transport model of main plasma, CDBM (current diffusive ballooning mode) model is applied.The diffusion coefficients by the CDBM model are given as follows [33]; where, c is the speed of light, ω pe is the electron plasma frequency, v A is the Alfvén velocity, q is the safety factor, s ≡ (r/q)(dq/dr) is the magnetic shear, α ≡ −(2µ 0 q 2 R/B 2 )(dp/dr) is the normalized pressure gradient.The fitting parameter C = 12.0, which is correction between the CDBM confinement time and the τ ITER 89 confinement scaling [34].The each fitting factor F(s, α), F κ [35] and F E × B [36] is defined as where is the E × B flow shearing rate and γ CDBM = |α| 0.5 (v A /qR)F(s, α) is the growth rate of CDBM [37].The radial electric field from the radial force balance is expressed as and the poloidal velocity V pol is evaluated with NCLASS.An additional thermal diffusion coefficient χ e = χ e + 0.1exp[−(ρ/0.075a) 2 ] is given for numerical stabilization, since the electron temperature can easily diverge due to the excessive radiation losses of impurities [38].The turbulent diffusion coefficient is set as D CDBM = χ CDBM /2 and turbulent convective velocity is set as V turb = − 0.5D turb r/a [39].

Impurity transport model
In the impurity transport calculations using the TASK/TI code (transport module for impurities), the following particle transport equations including atomic processes are solved for kvalent ionized impurity ions.
where, γ k is the ionization rate, α k is the recombination rate and S 0 is the impurity source.The impurity temperatures are assumed to be equivalent to the bulk ion temperature.The γ k , α k and radiation loss power of impurities are evaluated by using OPEN-ADAS database [40].Impurity densities at the boundary are given to be fixed value with the charge distribution determined by the coronal equilibrium.The time evolution of bulk plasmas and impurities is solved by linking TASK/TR and TASK/TI modules.The data of plasma density and temperature profiles are transferred from TR to TI, while the data of radiation loss power and effective charge number are transferred from TI to TR. Neoclassical transport coefficients for impurities are calculated by introducing geometrical factors (P A and P B ) into NCLASS in order to consider the poloidal asymmetry due to the toroidal rotation.The impurity flux is expressed in the following form [11,[41][42][43][44][45], with where the gradient length of each plasma quantity 1/L X = −∇X/X, and the aspect ratio ε = r/R.Absence of the poloidal asymmetry makes P A = 1 and P B = 0.The coefficients k i and C Z 0 are related to the collisionality, the neoclassical ion flow coefficient and the friction coefficient, respectively.For the case where the impurity is in the Pfirsch-Schlüter (PS) regime and the main ion is in the banana regime, − ( becomes close to −0.5, and a numerical formula is often used [46,47].More recently, extended formulas have been proposed with ion-electron collisional heat exchange for impurities, and C Z 0 is given as follows [48,49]: On the other hand, C Z 0 + k i is 0.33 takes the limit of ε ≪ 1, and the expression for k i is given by [50][51][52], which give the following expression including the effect of rotation [49], The variables f i , l i , g, α and µ ie in equations ( 17) and ( 18) are detailed in [48] and [49].Here, the trapped particle fraction f t is evaluated with Lin-Liu and Miller model [53]; Equation ( 14) assumes a convection coefficient K = 1 multiplied to the impurity density gradient and main ion density gradient terms.This assumption is known to be valid for the case where the impurity is in the PS regime and the main ion is in the banana regime.Therefore, rewriting equation ( 14) to equation (7) in [43].Using the temperature screening coefficient H, the neoclassical flux of the impurity is given as As shown in [43], the above analytical flux formula combined with the NCLASS code describes impurity transport including the poloidal asymmetry.Using the diffusion coefficient D NCLASS asym and convection velocity V NCLASS asym with the poloidal asymmetry (labeled with 'asym'), equation ( 20) can be simply expressed as where V NCLASS asym is defined as with the poloidally asymmetric convection velocity V NCLASS N,asym and V NCLASS T,asym .Here the poloidally asymmetric transport coefficients can be written as follows using the symmetric (labeled with 'sym') transport coefficients Therefore, the geometric factor that quantifies the poloidal asymmetry is , which increases the inward pinch driven by the main ion density gradient in the presence of strong asymmetry.Since from equation ( 14) and H is negative value except when the main ions are highly collisional regime, the screening effect driven by the main ion temperature gradient is reduced by the P B term in equation ( 25) [42].To evaluate the geometric factor P A and P B , a two-dimensional impurity distribution considering the poloidal asymmetry using flux surface averaged (FSA) density is given as follows [54]: Here, only the centrifugal effects due to toroidal rotation with NBI is considered, and the poloidally varying electrostatic potential caused by toroidal rotation-induced electrostatic potential, ion cyclotron resonance heating (ICRH) and electron cyclotron resonance heating (ECH) is not considered.
Assuming that the temperatures of the impurity species and bulk ions are comparable, the Mach number of impurities M Z is determined by M Z = √ m Z /m i M i , and Mach number of the main ion M i is determined by The toroidal rotation velocity V tor is a flux surface averaged quantity.The FSA value of the physical quantity X is calculated with where dl p is the line element along the poloidal magnetic field.In this simulation, the following turbulent diffusion coefficient, a combination of constant and parabolic ones, is used where, and the turbulent convection component of impurity is ignored.

Simulation of improved confinement plasmas
In this section, the results of transport simulations in the absence of impurities are presented using various JET plasma parameters to clarify the effects of improved confinement.Three cases are simulated: H-mode, hybrid scenario and advanced scenario plasmas.For the H-mode case, the plasma parameters of I p = 2.75 MA, B t = 2.5 T are used in [42].
Only NBI is applied as external heating and fuelling with P NBI = 17.5 MW at the radial position ρ NB 0 = 0.0 with the width ρ NB w = 0.6.The toroidal rotation velocity and pedestal profiles are given by the experimental profiles shown in [42].(JET#83351 at t = 12.37 s).The pedestal region is set for 0.9 ⩽ ρ ⩽ 1.0, and the radial profile is set with equation ( 5) for 0.8 ⩽ ρ ⩽ 1.0.The plasma density and temperature profiles are calculated by transport analysis for ρ < 0.8 without the pedestal region.In the hybrid scenario case, the plasma parameters of I p = 1.7 MA, B t = 2.0 T are used in [55].The NBI heating power profile is set to be P NBI = 16 MW with ρ NB 0 = 0.0, ρ NB w = 0.6, and the toroidal rotation velocity and the pedestal profiles are used in [55].In the advanced scenario case, the same plasma parameters as in the hybrid scenario cases are used for the NBI heating power, toroidal rotation velocity, and pedestal profiles, but a different value is used for the radial width ρ NB w = 0.8.For a current profile control, the additional LH heating with P LH = 2 MW is applied with ρ LH 0 = 0.4, ρ LH w = 0.1, and parallel refractive index N // = 2.1 [56].
Steady states for three cases are obtained.Figures 1 and 2 show the radial profiles of plasma quantities.As shown in figure 1, the safety factors of the hybrid scenario and advanced scenario plasmas have higher values in the entire radial region, compared to the H-mode plasma.At the plasma edge, q > 3 for the H-mode plasma, while q > 5 for the hybrid scenario and advanced scenario plasma.In the H-mode plasma, there is a region with q < 1, possibly causing sawtooth crashes, which are not considered in this simulation.In the advanced scenario plasma, the safety factor profile is flattened depending on the position of the LH heating, and negative magnetic shear regions are also formed.The bootstrap current fraction is f BS = 0.18 (H-mode), f BS = 0.21 (Hybrid scenario) and f BS = 0.26 (Advanced scenario).In the advanced scenario plasma, the fraction of ohmic current is less than 10% and the non-inductive current components are dominant.In the Hmode case, the plasma current density profile is monotonic, and is peaked at the magnetic axis where the plasma resistivity is small.In the hybrid and advanced cases, on the other hand, the plasma current is smaller than in the H-mode case, and the ratio of bootstrap current and externally-driven current are relatively large, resulting in a flattened or hollow current density profile.This gives a flattened or reversed q profile in the plasma confinement region, reducing the turbulent transport level to improve confinement.
In the H-mode case, the plasma density is higher than in the other conditions, and the profile is flattened for ρ < 0.8.The density gradient is relatively large in the plasma center due to the lower transport level and the stronger NBI particle source.The Mach number is M i = 0.26 at the plasma center, and a flat radial profile is given for ρ < 0.6.The H-factor and normalized-β value are H 98y2 = 0.87 and β N = 1.35, respectively.
In the hybrid scenario case, the plasma density is lower than in the H-mode case, but the plasma rotation is larger for 0.0 ⩽ ρ ⩽ 0.6, with M i = 0.36 at the plasma center and a maximum value of M i = 0.45 at ρ = 0.35.In the weak shear region inside ρ ∼ 0.4, turbulent transport is suppressed and the plasma gradient becomes stronger.Although a ITB is formed at ρ ∼ 0.4 in [42], the F E × B in equation (10) to be less than a few percent of the CDBM thermal diffusion coefficient, so that the effect of transport suppression due to toroidal rotation is small.A model for F E × B that can be adjusted based on experiments is proposed in [57].The H-factor and normalizedβ value are H 98y2 = 1.11 and β N = 1.97, respectively.
In the advanced scenario case, the plasma gradient increases steeply for ρ < 0.4 with the LH heating.The ITB formation results in larger temperature at the plasma center than that in the hybrid scenario plasma.This is due to the suppression of turbulent diffusion by the relaxation of magnetic shear s resulting from the LH current drive and by the increase of pressure gradient α resulting from the heating.The width of this increased region in the plasma gradient is approximately same with the width of the LH heating.The same toroidal rotation velocity is used, but the ion temperature is higher than in the hybrid case due to the formation of an ITB, giving a smaller Mach number M i = 0.30 at the plasma center.The Hfactor and normalized-β value are H 98y2 = 1.26 and β N = 2.34, respectively.

Validation of impurity transport model
We have introduced the impurity transport model described in section 2, and for the first step, the validity of the impurity transport simulation model is checked for the H-mode case in [42].For that purpose, calculations with fixed background plasma profiles are carried out.As shown in figure 3, the plasma density, temperature and toroidal rotation velocity profiles are used for two time slices at t = 9.32 s, when the W is not peaking in the center of the plasma (case A), and at t = 12.37 s, when the W is peaking in the plasma center region (case B).Here, it is assumed that the electron and ion densities are the same, and the electron and ion temperatures are the same.In case A, the plasma density is hollow in the plasma center, and has positive gradient.In case B, the plasma density is peaked in the plasma center, and has negative gradient.The toroidal rotation velocity profiles are flat in the plasma center with the Mach numbers of the main ion at the plasma center M i = 0.19 in case A and M i = 0.30 in case B. The W transport analyses using NEO and GKW show that the neoclassical transport is dominant in the plasma center, while the turbulent contribution is also significant in the plasma periphery [42].The turbulent diffusion component is nearly equal or larger than the neoclassical component for ρ > 0.45, and the turbulent convection component is larger and directing inward for ρ > 0.8.The turbulent diffusion coefficient of W is estimated from D W,GKW/ χ i,NEO in [42] to be set as in figures 4(a) and (b) by using equation (28), neglecting convective velocity.The neoclassical thermal diffusivity profile of the main ion obtained by NCLASS is used for the estimation assuming χ i,NEO = χ i,NCLASS .The density of W at the boundary is set as n w,edge = 2.5 × 10 14 m −3 in case A and n w,edge = 1.0 × 10 14 m −3 in case B.
Impurity profiles in the steady states in cases A and B are compared by using TASK/TI.The diffusion coefficient profiles are shown in figures 4(a) and (b).The neoclassical diffusion coefficient D NC is larger for ρ < 0.6 in case B because the higher density of main ion and W results in higher collisionality.In particular, the increase of D NC is significant for ρ < 0.15, where the W density profile is strongly peaking.As shown in figure 5, W does not accumulate in case A at the central region, whereas a strong accumulation occurs in case B. In case B, n W Z at the plasma center is about 1.5% of the electron density, so the quasi-neutrality condition is not significantly affected in fixed background plasma density.The densities at the low field side (LFS) in figures 6(c) and (d) are larger than the magnetic flux surface-averaged (FSA) W   Figure 7 shows that the normalized density gradient profiles of W are strongly depending on the magnitude and direction of the neoclassical convection velocities shown in figures 4(c) and (d).In both A and B cases, the neoclassical convective components V NCLASS T,asym driven by the main ion temperature gradient are positive in all the radial region.Furthermore, in case A, V NCLASS N,asym is positive for ρ < 0.6 due to the positive bulk density gradient, resulting in the absence of W density peaking at the plasma center.On the other hand, V NCLASS T,asym remains positive for ρ > 0.6, whereas V NCLASS N,asym becomes negative, and the difference of these terms gives negative V NCLASS asym for ρ > 0.7.As a result, inward convection is driven only in the outer region of the plasma, where a relatively small peak of the W density is formed.In case B, V NCLASS N,asym is negative in all the radial region, and V NCLASS asym is negative.In the outer region of the plasma, since the Mach numbers have similar value in both cases, the geometrical coefficients P A are similar, but at the plasma center, P A is increased to 40 in case B, while 10 in case A. Therefore, in case B, a sharp peak of the W density is formed at the plasma center.As shown in figures 6(a) and (b), the W density peak is shifted to the LFS.The poloidal asymmetry appears clearer in case A, when the central W density is lower.Note that in case B, the structure is also asymmetric.From figure 12 in [42], the reconstructed experimental values from SXR in case A are n W,LFS ∼ 5.6-9.2 × 10 15 m −3 at ρ ∼ 0.7 and n W,LFS ∼ 0.5-1.4× 10 15 m −3 at ρ ∼ 0.5, and are in good agreement with the results in figures 6(c) and (d).Similarly, the experimental values in case B are n W,LFS ∼ 0.3-1.1 × 10 16 m −3 at ρ ∼ 0.7 and n W,LFS ∼ 0.6-1.5 × 10 16 m −3 at ρ ∼ 0.15, and are also in good agreement.As shown in figure 7, the normalized W density gradients are similar to the results obtained with NEO and GKW in [42]; the difference around ρ = 0.3 in case A is due to the larger turbulent diffusion coefficient evaluated by GKW.

Tungsten transport in improved confinement plasma
Next, the impurity transport calculation is coupled with that of main plasma species.The cases with improve confinement plasmas in section 3 are analyzed to include the effects of impurity.In the H-mode case, the same W turbulent diffusion coefficient as in case B in section 4.1 are used.In the hybrid scenario case, the turbulent diffusion coefficient is estimated from [55] to have ρ W,turb 1 = 0.4, ρ W,turb 2 = 0.8, C W,turb 1 = 0.15, C W,turb 2 = 1.9, and C W,turb 3 = 1.5.In the advanced scenario plasma, ρ W,turb 1 = 0.6 is set from the previous conditions.The W density at the boundary n w,edge = 1.0 × 10 14 m −3 , which is same with the case in section 4.1, is used for all the cases.The toroidal rotation velocity profiles shown in figure 2(d) give central Mach numbers M Z ∼ 2.5 (H-mode), 3.5 (Hybrid), and 2.9 (Advanced).
In the H-mode case, peaked W density is obtained as show in figure 8.The normalized W density gradient in the center region of the H-mode case is similar to R/L nW ∼ 50 comparable to the case with fixed background which reproduces the result of the central accumulation.With the inclusion of W, the temperature gradient of the plasma is significantly reduced at ρ ∼ 0.1, due to stronger radiation losses with the impurity accumulation.Compared to figure 2(b), the electron temperature at ρ = 0.1 is reduced by about 4% and the normalized electron temperature gradient is reduced by about 40% due to the W accumulation.In this state, the dominant convection component is inward one driven by the bulk density gradient, as shown in figure 9(b).Thus, in this case of the W accumulation is large enough to cause a strong change in the electron temperature profile, and interactions occur between the main plasma and the impurity as follows: (1) increase radiation losses, (2) reduce main plasma temperature gradient and increase main plasma density gradient, (3) increase neoclassical inward flow, and (4) further W accumulation.A loop from (1) to ( 4) is followed until a steady state, resulting in enhanced and sustained W density peaking.Furthermore, the temperature screening effect is reduced by the term proportional to P B in equation (25) in the presence of a strong asymmetry in the W density.
In the hybrid scenario case, the W density also has a peak at the center of plasma.For ρ < 0.4, as shown in figure 10, the value of both V NCLASS T,asym and V NCLASS N,asym increase toward the center of plasma due to an increase in the main ion temperature gradient by the reduced turbulent diffusion and an increase in the main ion density gradient by the NBI particle source.The convective component driven by the main ion density gradient is stronger, resulting in an inward particle pinch for ρ < 0.35.The temperature screening effect is reduced in the center of the plasma, where the W density peaking is large, as in the H-mode case.Compared to the H-mode case, the W density is lower, whereas the W normalized density gradient is of similar value in the plasma center region.This difference is due to the strong inward neoclassical particle pinch driven in the H-mode case in the plasma periphery 0.85 < ρ < 1.0, which causes the inflow into the plasma confinement region with a significant amount of tungsten.In the advanced scenario case, as shown in figure 11, the particle pinch driven by the main ion density gradient is stronger due to the effects of reduced turbulent diffusion and the NBI particle source, as in the hybrid scenario case.Note that the neoclassical convection components are larger than in the hybrid scenario case, even though the Mach number is smaller, because the plasma transport is improved, and the plasma gradient is stronger than in the hybrid scenario case.The normalized W density gradient changes sharply at ρ ∼ 0.4, where the ITB is formed, and the W density peaks becomes wider in the radial direction than the hybrid scenario case.This result indicates the acceleration of W accumulation in the ITB region, as confirmed by previous experiments [14][15][16][17][18][19].

Parameter dependence of tungsten transport
In the previous section, it is shown that impurity accumulation occurs for improved confinement discharges.However, since impurity accumulation is not always observed in the experiments, we will investigate the conditions that lead to impurity accumulation from the dependence on the parameters.The first is the effect of rotation, which has a significant impact on impurity transport, and the second is the position of the ITB related to the performance of the transport barrier.

Toroidal rotation dependence
The dependence on Mach number of the main ion is investigated in the advanced scenario case to clarify the effect of toroidal rotation on the impurity profile formation.To avoid violation of the impurity trace limit condition α trace = n Z Z 2 /n i ≪ 1 and breakdown of calculation due to excessive radiation losses, the boundary condition of the impurity density is set to be smaller value n w,edge = 1.0 × 10 13 m −3 .Mach number M i0 at the plasma center is used here as a parameter, and the central toroidal velocity is set by V tor,0 = √ 2T i /m i M i0 .The profile of the toroidal velocity is given by using the function V tor =V tor,0 ( 1 − ρ 3 ) 2 to be similar to that of the hybrid case in figure 2(d).
In the small M i0 case as M i0 = 0.1, the W density profile flattens due to the turbulent diffusion without the enhacement by the toroidal rotation.The case with M i0 = 0.1 gives  P A = 3 at the plasma center, while the case with M i0 = 0.5 gives P A = 150, indicating formation of a strong asymmetric structure at the LFS. Figure 12 shows the dependence of the W density profile.The W density profiles in the high Mach number cases have the same kind of the structure as in the previous section.The W density profile flattens due to the turbulent diffusion, decreasing the neoclassical transport in the M i0 = 0.1 case.In addition, in the region near ρ ∼ 0.6, where outward W convection is driven, the enhancement by the toroidal rotation gives impurity ejection.The poloidal asymmetry of the W density distribution as in equation ( 26) scales to M 2 Z , and the W density is localized in the central region, so the value of P A is considerably larger at high Mach numbers.Figure 13 shows the dependence of the W normalized density gradient on the Mach number of W, plotting at ρ = 0.2, where the W density peaking is strong, ρ = 0.4, where it is close to the ITB-foot, and ρ = 0.6, where the outward particle pinch is driven.The value of the normalized W density gradient tends to decrease in accordance with decreased toroidal rotation velocity, especially at ρ = 0.2.Except for the higher Mach numbers, the values of V NCLASS N,asym and V NCLASS T,asym decrease with decreasing Mach number.At ρ = 0.6, the value of the W normalized gradient decreases from M Z ∼ 3-1 due to the reduction of enhanced temperature screening effect, but at ρ = 0.4, the value of the W normalized density gradient does not decrease due to a delicate balance of V NCLASS N,asym and V NCLASS T,asym .On the other hand, the normalized W density gradient is unchanged or even decreased at higher Mach numbers.This is due to the suppression of the temperature screening coefficient H/K via the M i dependence of k i as shown in [49].In [58], it is also theoretically shown that the W transport is suppressed in the high Mach number regime.

ITB position dependence
The dependence on ITB position is investigated in the advanced scenario case.The plasma parameters in section 4.2 are used, and the LH heating profile is changed with heating position ρ LH 0 = 0.4-0.7.The obtained plasma profiles show formation of the ITB, and the radial position of the ITB-foot is consistent with the LH heating position ρ LH 0 (figure 14).In the central region of the plasma, the ion temperature profile is peaking, but the electron temperature is flat due to the effect of radiation losses caused by W accumulation.The global plasma performance is improved to have H-factor H 98y2 = 1.39 and normalized-β value β N = 2.59 for ρ LH 0 = 0.6, compared to H 98y2 = 1.25 and β N = 2.39 for ρ LH 0 = 0.4.As shown in figure 15, the electron density at the plasma center is changed little with that at the ITB position, while the volume-averaged electron density increases from <n e > = 3.79 × 10 19 m −3 to <n e > = 3.85 × 10 19 m −3 for ρ LH 0 = 0.4-0.7.The W accumulation inside the ITB is suppressed in the cases of the ITB formed in the outer region.In the cases where the LH heating position ρ LH 0 are at the outer region, the particle source due to NBI is relatively small, so the outward V NCLASS T,asym driven by the main ion temperature gradient suppresses the inward particle pinch.On the other hand, in the cases where the   ρ LH 0 is close to the center of plasma, the inward V NCLASS N,asym causes a sharp increase in the W density gradient at the ITB region (figure 8).The concentration of W, c W = <n W >/<n e >, also decreases from 1.37 × 10 −5 (for ρ LH 0 = 0.4) to 0.89 × 10 −5 (for ρ LH 0 = 0.7).The W accumulation over the entire plasma can be suppressed by increasing the temperature screening effects.

Summary
Plasma transport simulations scheme including impurity transport have been developed in the integrated code TASK to predict the behavior of high-performance plasmas.The simulations are performed for the H-mode plasma, the hybrid scenario plasma, and the advanced scenario plasmas in a tokamak device.A neoclassical transport model including the effects of the toroidal rotation is applied to impurity transport analyses, because the large toroidal rotation is observed in mediumsized tokamaks.The calculated W profiles in H-mode plasmas show good agreements with experiments and first-principle simulations.In the plasma center, W accumulation can be enhanced due to the toroidal rotation effects, especially in the improved confinement mode where turbulent transport is suppressed.In addition, the accumulation of impurities becomes larger, because radiation losses cause an increase in the bulk plasma density gradient with a decrease in the bulk plasma temperature gradient.For smaller toroidal rotation velocity, the asymmetry of W is weakened, and the neoclassical transport is also reduced to suppress the central accumulation.Note that large toroidal rotation velocity can also enhance the neoclassical temperature screening effect, which gives a favorable impact on impurity ejection.In ITER, the central Mach number is predicted to be M i,0 ∼ 0.1 [49], so the effect of enhanced neoclassical transport of W due to toroidal rotation is expected to be smaller.On the other hand, the direction of particle convection is determined by the balance between the main ion density gradient and the temperature gradient, and impurity accumulation is possibly enhanced in the improved confinement mode plasma where the ITB is formed.The ITB plasma simulation in ITER predicts W accumulation in the plasma center [59], and it is necessary to control the plasma profiles to reduce the impurity accumulation.The present simulations show that the ITB formed at the outer region of the plasma can achieve higher plasma performance, including the suppression of impurity accumulation.
In this study, the transport simulation is carried out with focus on the neoclassical transport of W. However, it is known that turbulent transport and MHD activity have strong effects in the plasma periphery.In addition, it has been reported that minority heating with ICRH could enhance the neoclassical temperature screening effects [11], and a model is needed considering the background electrostatic field.In the future work, we will establish a simulation scheme including a turbulence model of impurity and simulate the entire plasma discharge including impurity control with RF heating.development research of the fusion DEMO reactor, and of RIAM in Kyushu University.Authors acknowledge discussions with Professor M. Honda.

Figure 1 .
Figure 1.Comparison of the radial profiles of safety factor q at the steady states of the H-mode, hybrid scenario and advanced scenario cases.

Figure 2 .
Figure 2. Comparison of the radial profiles of (a) plasma density, (b) electron temperature, (c) ion temperature, (d) toroidal rotation velocity, (e) Mach number of the main ion at the steady states of the H-mode, hybrid scenario and advanced scenario cases.The toroidal velocity profile in the advanced scenario case is the same as in the hybrid scenario case.

Figure 3 .
Figure 3. Referenced radial profiles of (a) plasma density, (b) normalized plasma density gradient, (c) ion temperature, (d) normalized ion temperature gradient, (e) toroidal velocity, (f ) Mach number of the main ion (solid line) and the W (dashed line).

Figure 4 .
Figure 4. Comparison of the radial W transport coefficients of turbulent and neoclassical diffusion in (a) case A and (b) case B, and neoclassical convective velocities in (c) case A and (d) case B.

Figure 5 .
Figure 5. Radial profiles of W density of the FSA value.

Figure 6 .
Figure 6.Calculated 2D plots of W density on the poloidal cross-section in (a) case A and (b) case B, and radial profiles of W density of LFS value on the midplane in (c) case A and (d) case B. The white lines in (a) and (b) indicate the equilibrium magnetic field configuration.The absolute values of the W density in (a) and (b) are significantly different, and the poloidal asymmetry is clearly visible in (a) since the central peaking is not occurring.

Figure 7 .
Figure 7. Radial profiles of normalized W density gradient of LFS value on the midplane.The dot represent the LFS value obtained from the NEO + GKW analyses in [42].

Figure 8 .
Figure 8.Comparison of the radial profiles of (a) W density and (b) normalized W density gradient with the H-mode, hybrid scenario and advanced cases.

Figure 9 .
Figure 9. Radial profiles of (a) W neoclassical convective velocities and (b) normalized main ion density and temperature gradients of the H-mode case.

Figure 10 .
Figure 10.Radial profiles of (a) W neoclassical convective velocities and (b) normalized main ion density and temperature gradients of the hybrid scenario case.

Figure 11 .
Figure 11.Radial profiles of (a) W neoclassical convective velocities and (b) normalized main ion density and temperature gradients of the advanced scenario case.

Figure 12 .
Figure 12.Comparison of the radial W density profiles, on the toroidal rotation velocities.

Figure 14 .
Figure 14.Change of the radial ion temperature profile depending on the LH heating position.

Figure 15 .
Figure 15.Dependence of the central W and electron density on the LH heating position.