The mixed gradient method [Zhong *et al.* Phys. Rev. Lett. **111**, 265001 (2013)] is adopted and effects of collisions and finite beta are included in the Weiland 9-equation fluid model. The particle flux and particle pinch, obtained using the Weiland anomalous transport fluid model, are compared with Tore Supra experimental results. Particle transport is also studied using predictive simulation data for an experimental advanced superconducting tokamak discharge in which neutral beam heating is utilized. The effects of collisions on particle transport are studied by turning collisions on and off in the Weiland model. It is found that the particle pinch region is related to the mode structure. The particle pinch region coincides with the region where the strong ballooning modes are present due to large gradients. The general properties of the fluid model are examined by finding regions where collisions can enhance the particle pinch.

## I. INTRODUCTION

It has been observed that central peaked density profiles can occur in tokamak discharges when edge fuelling is utilized.^{1–3} In Ref. 2, unresolved questions regarding particle transport are reviewed. Several studies concerning the possibility of inward particle transport (particle pinch) in tokamaks have been carried out an attempt to understand particle transport and particle pinches in present day experiments. Studies of the basis of particle pinches include an examination of neoclassical transport,^{4} an extension of the neoclassical study to include turbulence effects,^{5} pinches resulting from ion temperature gradient (ITG) drift waves,^{6,7} and the Weiland model describing strong particle pinch in the core due to electron trapping.^{8} The Weiland model provided the first description of a particle pinch that was strong enough to account for the density peaking in tokamak experiments with only edge fuelling.^{9} In particular, good agreement was obtained using the Weiland model to describe the evolution of the density profile in low collisionality discharges in the JET tokamak.^{10}

Exploration of efficient fuelling in the core of ITER tokamak plasmas is under investigation. A peaking of the density in ITER will depend upon the particle pinch due to drift waves.^{11,12} The occurrence of a particle pinch due to drift waves can improve the fusion energy gain factor, Q, for ITER. However, in ITER discharges, the effects of collisions may reduce the pinch effect.

A focus of the current paper is to understand and describe particle transport in tokamak discharges with medium collisionality. The effect of collisions on the development of a particle pinch was included in the studies described in Refs. 10 and 11 as well as those in Refs. 12–16. Dissipation due to collisions usually reduces the pinch effect. This is also true due to kinetic resonances. Thus, the fluid closure can significantly impact the results. This is clearly shown in the first quasilinear studies of ITG and Trapped electron (TE) modes,^{17–19} in which the goal was to obtain the particle pinch observed experimentally. These studies used expanded kinetic responses and concluded that the pinch they could obtain was much too weak to explain experimental results. This was later confirmed^{20,21} by using the unexpanded kinetic response.^{22}

Recently, an extended quasi-linear transport model, QuaLiKiz^{23} has been introduced. The QuaLiKiz model includes a frequency shift, which has been parameterized to the nonlinear frequency shift, obtained using a nonlinear gyrokinetic code. A particle pinch of the right order in magnitude is obtained using the QuaLiKiz model. It has been recently demonstrated^{24} that the nonlinear frequency shift in a fully nonlinear kinetic model motivates the fluid closure in the Weiland model. In the fluid limit it is, however, usually sufficient to keep only the linear part of the real part of the eigenfrequency. Zonal flows (not considered here) will enter as a nonlinear part added to the growth rate. Since it is kinetic resonances that remove the particle pinch in the standard quasilinear description, it is not surprising that QuaLiKiz model, with a frequency shift fitted to nonlinear gyrokinetics, yields a particle pinch of the right order of magnitude. The derivation in Ref. 23 involves fixed sources in configuration space. The fluid closure used in Refs. 7 and 23, as well as in this paper, is exact in the absence of sources in velocity space. This implies that the particle flux, calculated using the Weiland model, is valid for arbitrary curvature strength.

The advanced reactive toroidal fluid model, i.e., the Weiland model, both in its original local version (4-equation model) and in the more general version (9-equation model) is used to study particle transport in the Tore Supra and EAST^{1,3} (Experimental Advanced Superconducting Tokamak) tokamaks. The results obtained using the Weiland model are compared with the results presented in Ref. 23. In the local (4-equation model), ion continuity and ion energy equations, trapped electron density, and trapped electron energy equations are solved assuming Boltzmann free electrons, so it is an electrostatic version. In the more general treatment (9-equation model), parallel ion motion, impurity continuity and impurity energy equations, parallel momentum equation of free electrons, and Ampere's law are solved. Electromagnetic effects on all particle species are included and geometrical effects are taken into consideration.

The Weiland fluid model employed in this paper, and the kinetic model used in Ref. 25, yield results that are surprisingly similar. The similarity in results is due to the “frequency width” in the QualiKiz model where it effectively removes Landau damping, thus restoring the particle pinch. However, there are some significant differences as well.

One difference is that with the 9-equation model, unstable electron modes are found for all gradients considered. This difference may be due to the presence of the interchange drive in the 9-equation model, in which electromagnetic effects are included. These electromagnetic effects are ignored in Ref. 25. However, the strongest electron mode in the Weiland model appears in the region where QualiKiz also has a TE mode. A second difference is that the experimental pinch data from Tore Supra follow the 4-equation electrostatic, collisionless separation between the stability regions more closely than obtained in Ref. 25 as shown in Fig. 6 of Ref. 25. Another important difference is that with the Weiland model, the pinch region is found to be coupled to the mode structure such that the pinch region coincides with the region of strong ballooning (well localized eigenfunction) when the gradients are increased.

## II. FORMULATION

As noted in Ref. 25, the sign of the flux depends on a linear combination of temperature and density gradients. The convective terms due to the magnetic field gradient are present in fluid models considered in this paper, but these terms will be kept fixed in order to focus on density and temperature gradient effects. Usually, a competition between convection and compression is observed. However, for the parameters used here, the flux will have one part proportional to $R/Ln$ and the other proportional to $R/LT$, where *R* is major radius, and $Ln=\u2212(dlnn/dr)\u22121$ and $LT=\u2212(dlnT/dr)\u22121$ are scale lengths of density and temperature inhomogeneity. The assumption is made that the ion and electron temperature length scales vary together, that is, $LT=LTi=LTe$.

It is convenient to introduce a linear combination of density and temperature gradients of the form

Note that the linear combination used in Eq. (1) differs from the linear combination used in Ref. 25. This combined gradient results in a straight line in the $(R/Ln,R/LT)$ plane. It is found that employing the advanced reactive model^{26} in the local 4-equation limit, the result obtained for the regions of stability and instability of the ion temperature gradient and trapped electron modes, shown in Fig. 1, is very similar to the results presented in Fig. 1 of Ref. 27.

The results presented in Fig. 1 in this paper are obtained taking $Te/Ti=1$ and $FT=0.356$, where FT is the fraction of trapped electrons. Consequently, $\zeta =R/LT+0.68R/Ln$ and the critical gradient is *ζ* = 10.5. The inclusion of Fig. 1 provides an introduction and overview of the problem.

The stability regions plotted in Fig. 1 are found to be similar to the stability regions shown in Fig. 6(a) in Ref. 25 (included here as Fig. 3(c) for comparison). A particle pinch is found in the entire plane except near marginal stability of the ITG mode when the 4-equation model is utilized. Actually, for the parameters used here also the trapped electron mode gives a particle pinch. The collisionless particle flux given in Eq. (2) below was obtained using the 4-equation model presented in Eq. (6.158) of Ref. 21

where $\u03f5n=2Ln/R,\u2009\eta e=Ln/LT$, and $\omega \u0302$ is the real frequency normalized by magnetic drift frequency. The *η _{e}* component in the second term on the right side of Eq. (2) gives a particle pinch for the ITG modes $\omega \u0302r<0$. The

*η*component in the last term on the right side of Eq. (2) can provide a pinch as well, even for the trapped electron mode (positive), where the part in the second term gives an outward flux, the last term dominates already for $\u03f5n=0.7$ then giving a net pinch for the trapped electron mode. Thus, the conventional wisdom that ITG gives a pinch and the TE mode yields the outward flux is an oversimplification.

_{e}Note that Eq. (2) is from the original derivation^{9} so the fact that the TE mode can give a particle pinch has been available since 1989. Since the pinch due to the trapped electron mode occurs for substantial values of *ϵ*_{n}, it is not leading to a very peaked density profile. It does, however, contribute to the general picture that inhomogeneity length scales in different fields tend to become of the same order and certainly may contribute dynamically in the initial phase when a flat density profile becomes peaked. The dynamic phase may also enter in transient transport experiments.

The equation of particle flux and effective diffusivity *D*_{eff} can be written as by employing Ficks law

The effective diffusivity can also be written in matrix component form

where *D _{n}* is the diagonal part associated directly with density peaking;

*D*is the off diagonal part;

_{e}*V*is the convective part; and

*GK*is the result of a geometry effect. The variables DDIAG and VCONV are defined as $DDIAG=Dn$ and $VCONV=Deff\u2212DDIAG$. Note that all diffusion and convection coefficients here depend on the linear eigenvalue

*ω*and, as a consequence of that dependence, on all gradients. Therefore, it is not possible to obtain a complete separation of all parts analytically. In fact, the picture would change if

*ω*is expressed through the dispersion relation. Nevertheless, the different coefficients obtained in Eq. (4) appear in a natural way in the derivation and thus they are displayed separately in Sec. III.

## III. RESULTS OBTAINED APPLYING TORE SUPRA CONDITIONS

The stability region and the particle flux in the $R/LT,R/Ln$ plane, calculated utilizing the 9-equation model, are shown in Figs. 2(a) and 2(b), respectively. Although the momentum transport is not studied here, the 9-equation model includes parallel ion motion with general degree of ballooning, magnetic-q and magnetic shear effects, electromagnetic effects, and the effects of electron-ion collisions on trapped electrons.

If the line with $R/Ln=2.06$ is followed and the effective diffusivity computed with $\eta e=\eta i=\eta $, one obtains the results for *D*_{eff} as function of *ζ* that are shown in Fig. 3(a). Note, a particle pinch is found between $\zeta =8.8$ and $\zeta =9.6$. As discussed below, the pinch occurs in a region where there is strong ballooning in the solution of the linear eigenvalue problem. A Gaussian eigenfunction is assumed in the model of the form

The degree of ballooning is given by the exponent *α* which is particularly large in the pinch region as seen in Fig. 3(b).

A fairly rapid variation of the exponent around *ζ* = 9 is found where the pinch gets strong. The increase in the exponent is due to a feedback mechanism, when the growth rate increases the localization increases, which makes the growth rate to increase more. For the parameters used here, the ITG gives a pinch, while the TE mode gives outward transport. In the region of weak ballooning, the curvature drive of the ITG mode with average space dependence $<cos\u2009\theta +s\theta \u2009sin\u2009\theta >$ tends to be averaged out over the mode profile so the ITG gets weak, while the TE mode remains strong and dominates. In the strong ballooning regime, the ITG gets much stronger, dominates, and creates a particle pinch.

The reason why the localization is then reduced in Fig. 3(b) is that the ITG mode changes the sign of the real part of the frequency and becomes an electron mode. Consequently, this electron mode dominates for *ζ* > 10 which happens to be the region in Fig. 1 where there are both ITG and TE modes. A similarity of the 9-equation model, for *ζ* > 10, with the 4-equation model in Fig. 1 is that the local approximation is valid. Since the main ITG mode turns into an electron mode in this region, the particle flux is directed strongly outward and this is where the asymptote in Fig. 3(a) is found.

These results presented here are found to be consistent with the results shown in Fig. 6 of Ref. 25 (included as Fig. 3(c) below). Note that the 9-equation model applied in this paper is electromagnetic, while the model used to obtain the results presented in Ref. 25 is an electrostatic model. This particle pinch is found where the experimental circle-dashed dots are located in Fig. 6 of Ref. 25. The asymptote (large outward flux) appears at $\zeta =10.1$. Thus, the Tore Supra experimental particle pinch is located in the region where the 9-equation model indicates the occurrence of a pinch.

It is interesting to compare the results obtained using the 9-equation model with the results obtained in Ref. 25. The mixed gradient $\zeta =R/LT+0.68R/Ln$ is used in the 9-equation model, since *ζ* is constant along the tangent to the boundary between ITG and ITG + TE unstable regions shown in Fig. 1. This constant value is 10.5. However, in Ref. 25, the authors used a linear combination which is constant on the corresponding boundary obtained with QuaLiKiz with $\zeta c=R/LT+4R/Ln$ and the tangent is $\zeta c=22$.

The high outward flux region in Fig. 6 of Ref. 25 is located at $\zeta c>21.7$, while for the calculations presented here it is at $\zeta c>19.7$ for $R/Ln=2.8$. Thus, the results are quite close. The largest difference appears for smaller $R/Ln$ where QuaLiKiz requires $R/LT\u224812$ for outward flux, while the Weiland 9-equation model only requires $R/LT=8$.

The differences between the 4-equation model and the 9-equation model are in the treatment of collisions, electromagnetic effects, and as well as effects associated with solving the eigenvalue problem, in particular, the degree of localization of the eigenfunction. The nonlocalized eigenfunction is associated with weak ballooning; the well localized eigenfunction, with strong ballooning; and the limit of very strongly localized eigenfunction, the local limit. It is interesting to include the 4-equation model for several reasons. First, the stability areas in Fig. 1 and Fig. 6 of Ref. 25 are quite similar. The stability regions obtained with the 4-equation model are also in close agreement with those shown in Fig. 1 of Ref. 27. Second, it turns out that the maximum pinch region, as shown in Fig. 6 from Ref. 25 (circle-dashed dots), is much closer to the onset of the TE mode in Fig. 1 than in Ref. 25. Also the slope of this boundary is much closer to the slope of the strong pinch region in the results presented in this paper. The authors of Ref. 25 indicate that “the density and temperature gradients in steady state should be very close to the TEM boundary.” This is true for the Weiland 4- and 9-equation models, but not for the QuaLiKiz model, at least not for $R/LT=2.2$.

It is found that collisions make the transport more outward both for small and large *ζ* as shown in Fig. 4. However, there may be regions with the opposite trend for intermediate *ζ*. The effective collision frequency of experimental data for trapped electrons is $\nu \u2009exp\u2009\u22480.4$. It is seen in Fig. 4 that the sensitivity to electromagnetic effects may be stronger than that on collisions. As noted above, there is also a pinch due to the TE mode in the collisionless case and that is contributing strongly to the pinch in Fig. 4(b). The general picture is that there is a dominant electron mode everywhere except near the pinch region in the electromagnetic case. The transition from ITG to TE mode occurs through a smooth change in the sign of *ω* from negative to positive when we pass into the region of local modes (*ζ* = 10). However, in the electrostatic case, the ITG mode dominates in the whole plane below *ζ* = 10.

The diagonal and convective fluxes, defined above in Eq. (4) and the paragraph following in Eq. (4), are shown in Fig. 5 as a function of the mixed gradient variable *ζ* for the Tore Supra discharge corresponding to Fig. 3(a). DDIAG is the diagonal component of the diffusivity associated directly with density peaking.

## IV. EAST TOKAMAK PINCH RESULTS

A simulation of an EAST tokamak discharge heated with neutral beam injection yields parameters that are then used in the study of pinch effects in the EAST tokamak. The resulting parameters are: $Te/Ti=1.8,\u2009FT(Trapped\u2009Particle Fraction)=0.36,\u2009Te=2.8\u2009keV,\u2009n=3\xd71019\u2009m\u22123$, *q* = 3, *s* = 1, $\nu \u0302=0.4$, and $Btor=1.5T$. Note, the collisionality is comparable with that for Tore Supra. The model employed in carrying out the pinch study is the 9-equation model which includes finite *β* effects. It is found that the mixed linear gradient threshold is

By following a line with constant $Ln\u2009(R/Ln=4.11)$, the diffusivity, as a function of *ζ*, is obtained and shown in Fig. 6(a). The negative values of *D*_{eff} in the region between *ζ* = 4.6 and *ζ* = 5.1 indicate a density pinch.

In order to illustrate the effect associated with the eigenfunction, the diffusivity is obtained with *q* = 25. When *q* is increased to 25, corresponding to strong ballooning, the result obtained is shown in Fig. 6(b). In this case, the linear combination threshold is

The plot of *D*_{eff} in Fig. 6(b) indicates a stronger pinch effect over a broader range of *ζ*. For Tore Supra, it was found that the pinch occurs in the strong ballooning limit. Thus, it is not surprising that a similar result is obtained for the East discharge with *q* = 25. In Fig. 7, the analytical diagonal and convective transport coefficients are shown for the EAST conditions used in obtaining the result shown in Fig. 6(a).

Note that the asymptote in EAST result appears for lower *ζ*, *ζ* = 5.2 than in Tore Supra, *ζ* = 10.2. This is a consequence of the impact of the different aspect ratios of the two tokamaks on the value of $R/Ln$. Although, for Tore Supra, the line *ζ* = 10.5 in Fig. 1 is obtained applying the two dimensional 4-equation model that has only toroidal drive, it is noted that the parallel ion motion becomes subdominant for large linear growth rate. This is the local limit. The increased outward flux at the gradient where the TE mode enters in the 4-equation model may be seen as a consequence. However, the TE mode also results in a pinch above *ζ* = 10.5 in the application of the 4-equation model. The collisions on trapped electrons and electromagnetic effects result in outward transport under these conditions.

The result when the 9-equation model is used is that for small gradients, the ITG mode is in the weak ballooning region where the driving curvature term is almost averaged out. In this region the TE mode dominates because it is evaluated for deep trapping which then implies that the driving curvature term is evaluated on the outside of the tokamak. For larger *η*, the ITG mode comes into the strong ballooning region and takes over from the TE mode. As a result the pinch occurs. Then, at even higher *η* the electron mode takes over both in the Tore Supra and EAST cases.

## V. CONCLUSIONS

It is shown that collisions and finite beta play important roles with regard to the direction of the particle flux, both in Tore Supra and in EAST discharges. The fluxes in the simple 4-equation model, giving the mode boundaries in a $R/LT,R/Ln$ diagram are compared and found similar to the results shown Fig. 6(a) in Ref. 25. The fluxes obtained using the 9-equation model are nearly the same as those Fig. 1 of Ref. 22. It is observed that collision and finite beta effects destabilize the electron modes in almost the entire $R/LT,R/Ln$ plane. These effects also change the direction of the particle flux from inward to outward for mixed gradients above the collisionless separation line *ζ* = 10.5. Another feature of the 9-equation model is that it includes parallel ion motion with an arbitrary degree of ballooning. This leads to dominance of the trapped electron mode in a major part of the domain. However, an increased growth rate turns the ITG mode into the strong ballooning limit which makes it dominant. This leads to a net particle pinch quite close to but below the collisionless separation line *ζ* = 10.5. The effect of strong ballooning was shown for using parameters obtained in a simulation of neutral beam heated EAST discharge. The 4-equation model is two dimensional and includes only the toroidal drive of ITG. Since the toroidal drive dominates for large *η*, the 4-equation model may be more relevant for large *ζ*. It is then only necessary to add collisions and electromagnetic effects^{28–30} to recover outgoing flux as obtained by applying the 9-equation model. In this work, only the total flux has been compared with experiment but diagonal and convective fluxes have been discussed according to their analytical definitions. With utilization of this definition, there is no convective flux reversal resulting in a pinch found in the results presented in this paper. The reversal of the total flux is due to the balance between diagonal and convective fluxes as a whole. Both diagonal and convective fluxes depend on the linear eigenvalue and that depend on all gradients in the system.

## ACKNOWLEDGMENTS

We are grateful to W. L. Zhong, S. D. Song, and X. L. Zou for discussions on Ref. 25. This work was partly supported by “One-Three-Five” Strategic Planning Project of Chinese Academy of Science No. Y25FZ10051 and partly supported by Chinese Academy of Sciences visiting professorship for senior international scientists No. 2012T1J0026. Two coauthors, T. Rafiq and A. H. Kritz, acknowledge support by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Contract No. DE-FG02-92-ER54141.