# Integrated modelling of neon impact on JET H-mode core plasmas

M. Marin, J. Citrin, C. Giroud, C. Bourdelle, Yann Camenen, L. Garzotti, A.

Ho, M. Sertoli, Jet Contributors

► To cite this version:

M. Marin, J. Citrin, C. Giroud, C. Bourdelle, Yann Camenen, et al. Integrated modelling of neon impact on JET H-mode core plasmas. Nuclear Fusion, 2022, 63 (1), pp.016019. 10.1088/1741-4326/aca469. hal-04234737

HAL Id: hal-04234737

https://hal.science/hal-04234737v1

Submitted on 10 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

# Integrated modelling of neon impact on JET H-mode core plasmas

To cite this article: M. Marin et al 2023 *Nucl. Fusion* **63** 016019

View the [article online](#) for updates and enhancements.

You may also like

- *Flux-driven integrated modelling of main ion pressure and trace tungsten transport in ASDEX Upgrade*
O. Linder, J. Citrin, G.M.D. Hogeweij et al.

- *Core turbulent transport in tokamak plasmas: bridging theory and experiment with QuaLiKiz*
C Bourdelle, J Citrin, B Baiocchi et al.

- *Multiple-isotope pellet cycles captured by turbulent transport modelling in the JET tokamak*
M. Marin, J. Citrin, L. Garzotti et al.
# Integrated modelling of neon impact on JET H-mode core plasmas

M. Marin<fnref n="1" />,*, J. Citrin<fnref n="2" /> ID, C. Giroud<fnref n="3" />, C. Bourdelle<fnref n="4" /> ID, Y. Camenen<fnref n="5" />, L. Garzotti<fnref n="3" /> ID, A. Ho<fnref n="2" /> ID,
M. Sertoli<fnref n="6" /> and JET Contributors<fnref n="1" />

E-mail: michele.marin@epfl.ch

Received 10 August 2022, revised 18 November 2022

Accepted for publication 21 November 2022

Published 12 December 2022

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

Nuclear fusion reactor plasmas will need to exhaust a significant proportion of energy flux through radiative processes, to enable acceptable divertor loads. This can be obtained by line radiation from impurities, injected from the plasma edge. There are however limitations on the sustainable impurity content, since radiation from the core can lead to a deleterious electron heat sink. Moreover, dilution of the main ions reduces the available fuel. Simultaneously, impurities have an impact on the turbulent transport, both by dilution and by changes in the effective charge. Recent experiments at JET point towards an improvement in plasma confinement in neon seeded discharges with respect to purer equivalent plasmas. In this paper the impact of the impurities on the confinement is studied, isolating various effects.

First-principle-based integrated modelling with the QuaLiKiz quasilinear turbulent transport model explains the improvement by a combination of higher pedestal temperature, increased rotation shear, and impurity-induced microturbulence stabilization. These results are optimistic with respect to the maximum impurity levels allowed in ITER and future reactors. Comparison between QuaLiKiz and higher fidelity gyrokinetics has exposed issues with QuaLiKiz impurity peaking predictions with rotation.

Keywords: integrated modelling, seeding, turbulence, first principles

(Some figures may appear in colour only in the online journal)

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹ SPC, Swiss Plasma Center, Ecole Polytechnique Federale de Lausanne, Station 13, Lausanne 1015, Switzerland</footnote>

<footnote n="2">² DIFFER—Dutch Institute for Fundamental Energy Research, Eindhoven, Netherlands</footnote>

<footnote n="3">³ CCFE, Culham Science Centre, Abingdon OX14 3DB, United Kingdom of Great Britain and Northern Ireland</footnote>

<footnote n="4">⁴ CEA, IRFM, F-13128 Saint-Paul-lez-Durance, France</footnote>

<footnote n="5">⁵ CNRS, Aix-Marseille Univ., PIIM UMR, Marseille, France</footnote>

<footnote n="6">⁶ Tokamak Energy Ltd, Culham Science Centre, Abingdon OX14 3DB, United Kingdom of Great Britain and Northern Ireland</footnote>

<footnote>i See Mailloux et al 2022 (https://doi.org/10.1088/1741-4326/ac47b4) for JET Contributors.</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

Heat fluxes to the divertor in ITER and in future reactors will need to be mitigated. The peak heat load that can be managed by the actively cooled ITER divertor is $q_⊥ \sim 10 \text{ MW m}^{-2}$ [1]. During the International Thermonuclear Experimental Reactor (ITER) full-power fusion operation, the power from the confined plasma to the scrape off layer (SOL) is expected to be ~100 MW. The peak heat flux at the divertor depends on the SOL power fall-off length $λ_q$, projected to be as low as ~1 mm [2]. $q_⊥$ will then depend on the exact geometry of the equilibrium, but with such low $λ_q$ the operational limit is expected to be exceeded. One of the strategies to reduce the heat load on the divertor is through impurity radiation. Even though there is a large uncertainty in the projections for the value of $λ_q$ in ITER, a minimum of ~30 MW of light impurity radiation from the divertor region is foreseen to be required.

In the Japanese DEMOnstration power plant (DEMO) a $P_{rad}/P_{tot}$ of 0.8 is envisioned [3], while for the European DEMO a value as large as 0.9 might be necessary [4]. In present day experiments, significant radiation fractions have been obtained with impurity seeding while maintaining good confinement and compatibility with edge localized modes (ELMs)-free regimes. For example, in the Axially Symmetric Divertor Experiment Upgrade (ASDEX Upgrade), significant edge radiation with argon seeding was achieved concomitant to ELM suppression with Resonant Magnetic Perturbations [5].

Low Z impurities have the advantage of radiating mostly in the 10–100 eV range, while almost exclusively Bremsstrahlung is emitted from the core, where they are fully stripped. In the Joint European Torus (JET), the reduction of the heat loads on the divertor was demonstrated in a variety of divertor configuration both in carbon wall [6] and in ITER like wall (ILW) [7]. Nitrogen is one of the most used and studied element for impurity seeding [8–10], with effects on the heat flux and, importantly, on the turbulent transport.

Improvement of confinement, and of the pedestal in particular, has been observed both in high triangularity, ELMy H-modes at JET [11] and ASDEX-U [12]. Multiple effects could contribute to this change in the pedestal, but a quantitative model to predict the profiles in this region is not available.

Considering the effect of impurities in the core, there are further disadvantages, most notably dilution and core radiation. Dilution decreases the amount of fuel in the plasma, while radiation decreases the temperature and the heat flux at the separatrix. Tungsten influx, for example, can deteriorate the pedestal and alter the frequency of the ELMs. ELMs regulate the impurity content in the core by flushing out particles [13], so a decrease in temperature can be followed by further impurity accumulation [14]. Depending on divertor conditions, this might create a feedback loop and cause a disruption.

The use of Nitrogen in ITER might prove to be problematic, mainly due to the formation of ammonia [15], which would limit the achievable duty cycle [1]. An additional candidate as seeded impurity is neon. In JET carbon and ITER-like walls neon is an effective radiator [16, 17], but the radiation is localized further in the pedestal region with respect to Nitrogen,

which often results in loss of pedestal performance. However, in ITER neon radiation is predicted to be more compressed in the divertor region, which might modify this picture.

In the Tokamak Experiment for Technology Oriented Research and DIII-D, improvement in confinement was observed with Neon seeding in L-mode plasmas [18–20]. Degradation of the pedestal was instead observed in high triangularity ($\delta \sim 0.4$) experiments, resulting in an overall decrease of the confinement [6]. Similarly, in ILW, $\delta = 0.4$, relatively low power ($P_{NBI} \sim 20 \text{ MW}$) discharges, neon was observed to be detrimental. At higher power, the pedestal degradation was mostly balanced by an increase in core performance [21]. This lower core transport was observed at $\delta = 0.2$ and $Z_{eff} \sim 1.8$.

More recent scenarios for the C38 campaign at JET, conversely, retrieved the enhancement in confinement observed in Nitrogen. The aim of the experiments was to compare Ne and N as radiators in JET as a benchmark for ITER. The discharges had high power, $P_{NBI} = 25 \text{ MW}$, triangularity, $\delta = 0.4$ and reached $Z_{eff} \sim 2.5$ The improvement was significant and found to correlate with the neon content and to require exceeding an input power threshold. Interpretive analysis showed improvement of both pedestal and core transport.

Impurity transport has been extensively studied [22]. Nitrogen seeding has an impact on the core turbulent transport, as was observed and calculated for L-mode plasmas at JET [23]. The improvement was found to originate from the change in $Z_{eff}$, which in turn caused a change in the ratio between the magnetic shear and the safety factor $s/q$ and a different $T_e/T_i$, where $T_e$ and $T_i$ are the electron and ion temperatures. In neon seeded, carbon wall discharges, the increase in $Z_{eff}$ was shown to lead to higher core $T_e$ and $T_i$ [24]. Impurity density peaking is also linked to turbulent stabilization [25], but high quality impurity profiles are not easily available and require dedicated experimental analysis, despite recent progress [26]. At the same puffing levels, the core transport reduction with neon was observed to be larger than with nitrogen. The larger $Z_{eff}$ of the neon plasmas at constant dilution is a candidate in explaining this effect.

Previous simulations for ITER [27, 28] predicted only a minor decrease in the alpha power following neon seeding, due to dilution being partially balanced by the higher ion temperature. Both simulations calculated or assumed a similar pedestal in the seeded case. Understanding the precise origin of the improvement in confinement is therefore important to increase confidence in extrapolability and justify further pedestal analysis. Despite recent progress [29], the study of the changes in the pedestal is very challenging, due to the lack of a quantitatively accurate and sufficiently fast pedestal model. Conversely, the tools to study the differences in the core transport are available. Integrated modelling with reduced-order turbulent transport models such as QuaLiKiz, which can handle an unlimited number of ion species, [30] is an ideal framework to tackle this problem and allows to quantitatively reproduce the experiments.

The paper is structured as follow. Section 2.1 presents the two JET discharges that were compared, with and without high neon seeding. In section 2.2 the integrated modelling settings are presented. Section 2.3 explains the results of the modelling
and breaks up the impact of the various differences. Section 3 explores the sensitivities of the modelling, aiming at isolating the impact of the impurities on the kinetic profiles. Section 4.1 shows the QuaLiKiz standalone predictions and sensitivities and section 4.2 compares the QuaLiKiz predictions with results from the higher fidelity gyrokinetic code GENE [31]. Conclusions are discussed in section 5.

## 2. Unseeded and neon—seeded comparison

### 2.1. Experimental discharges

Neon seeding was studied during the recent C38 experimental campaign at JET. A scan was performed, from discharge #96132 to #96139, with decreasing neon puffing. The analysis in this paper will consider the two extremes, unseeded discharge #96139 and high seeding #96133, both in deuterium. Average $Z_{eff}$ was ~ 1.4 and ~ 1.8 respectively. The details are presented in table 1. The experiments were performed at high triangularity and ITER-relevant, vertical target divertor configuration. The electron density $n_e$, $T_e$, $T_i$ and the toroidal rotation $\omega$ were extracted and fitted using the Gaussian process regression tools available at JET [32]. The experimental data were fitted on the 12.5–13 s and 14–14.5 s time intervals for the unseeded and the seeded cases respectively. The equilibria were reconstructed by pressure-constrained EFIT++ [33, 34]. The whole equilibria were radially shifted by 0.63 cm and 0.0 cm for #96133 to #96139 respectively. This was done to ensure an electron temperature measured by the High Resolution Thomson Scattering (HRTS) of 100 eV at the separatrix, as predicted through power balance by the two point model illustrated in [35], and as done in previous works [36, 37]. The density profiles were scaled by 1.0 and 1.007 in order to be consistent with line integrated interferometry measurements. As can be seen in figures 1 and 2, the measurements from Laser Imaging Detection and Ranging (LIDAR) have large errorbars for this discharge. The laser energy in these pulses was around half the optimum, which led to poor performance of the diagnostic. While comparing to the simulated results, HRTS should be trusted over LIDAR. The latter is still kept since it provides information about the shape of the profiles for $\rho < 0.25$, where HRTS is not available. $\rho$ is the normalised toroidal flux coordinate $\rho_{tor} = (\frac{\psi_{tor}}{\psi_{tor,LCFS}})^{\frac{1}{2}}$.

For the ion temperature and rotation measurements, the neon seeding allows for high throughput of the neon charge exchange (CX). The neon channels looking at 524.9 nm, the NeX line, were chosen to extract the core data. Due to significant scatter between the CX channels looking at different wavelengths, the other measurements were removed from the fitting and are not showed in the comparisons. Sawteeth were present with similar magnitude in both discharges, with an inversion radius and frequency of $\rho \sim 0.25, f \sim 2$Hz for #96133 and $\rho \sim 0.3, f \sim 4$Hz for #96139.

### 2.2. Modelling settings

The integrated modelling is performed using the JETTO transport solver, within the JINTRAC [38, 39] suite of codes.

While JINTRAC has the capability of modelling both edge and core, only the latter was evolved. Within JINTRAC this is handled by JETTO, which is a 1.5 dimensions transport solver. The Faraday's equation for the time evolution of the plasma current, the electron and total-ion energy equations and the individual ion mass-continuity and total ion toroidal momentum equations are solved. A number of modules calculating sources, sinks and transport are called with various frequencies, maximizing the speed and maintaining the accuracy. The relevant modules are mentioned below. In this work, the current profile $j$, $T_e$, $T_i$, and $n_e$, are predicted self consistently. The experimentally measured rotation profile is prescribed. The impurities densities $n_{imp}$ are either imposed from the experimental measurement or self consistently evolved, depending on the set of simulations. Four impurities are included: beryllium, nickel, tungsten and neon (when puffed). At the beginning of the predictive simulations and through the entirety of the simulations where impurities are imposed, the densities are set to match

* The line-of-sight integrated measurement of $Z_{eff}$

* The soft x-ray emission

* The observed poloidal asymmetry of the soft x-ray emission

* The line-of-sight integrals of the total radiation as measured by bolometry

The method presented in [40] was applied to infer the impurity profiles. The simulations were run for two seconds of discharge time, which corresponds to ~10 energy confinement times, sufficient to reach a stationary state.

The equilibrium is evolved using ESCO [39], with the resistivity calculated by NCLASS [41] and the boundary conditions prescribed from EFIT++. The neutral beam injection (NBI) heating is modelled using PENCIL [42], while PION [43] is used for the ion cyclotron resonance heating (ICRH) and for a more accurate approximation of the fast ion slowing down time. SANCO [44] is utilized to calculate the charged states of impurity and the radiation, assuming coronal equilibrium as initial condition. No attempt is made to model the pedestal and internal boundary conditions are imposed at $\rho = 0.85$.

The neutral source is modelled by FRANTIC [45] and the particle flux crossing the separatrix is set to $10^{21} \text{ el s}^{-1}$. The experimental target for the D particle flux was $2.6 \times 10^{22} \text{ el s}^{-1}$ and a penetration efficiency on the order of 10% is normally assumed for the modelling. Due to the relatively high power, a lower value was chosen. There is no consistent methodology to measure the penetration efficiency, but since only a fraction of particles are deposited inside $\rho = 0.85$ a low sensitivity can be expected. Dedicated simulations confirmed an extremely low impact of the profiles to the choice made for the neutrals source. The neutron production rate, mostly due to beam-thermal reactions is calculated within PION. The turbulent transport for main species and impurities is predicted by QuaLiKiz, while NCLASS is used for the neoclassical transport. Simply summing neoclassical and turbulent transport for impurities is justified for light impurities and relatively low collisionality [46].
<table><caption>**Table 1.** Key parameters of JET shot #96139 and #96133. $I_p$ and $B$ are the total plasma current and the on-axis magnetic field, $P$ represent the injected heating power. $\beta_N$ is the normalized plasma beta, defined as $2\mu_0^2 a_{min} * W_{tot}/(I_p * V * B_{geo})$, with $a_{min} = 0.5 * (R_{out,LCFS} - R_{in,LCFS})$, $B_{geo}$ the vacuum toroidal field at the geometric plasma centre and $V$ the volume.</caption><thead><tr><th>Shot #</th><th>$I_p$ (MA)</th><th>$B$ (T)</th><th>$Z_{eff}$</th><th>$P_{NBI}$ (MW)</th><th>$P_{ICRH}$ (MW)</th><th>$\Phi_D$ (10<sup>22</sup> at s<sup>-1</sup>)</th><th>$\Phi_{Neon}$ (10<sup>22</sup> at s<sup>-1</sup>)</th><th>$\beta_N$</th></tr></thead><tbody><tr><td>96139</td><td>1.4</td><td>1.7</td><td>1.4</td><td>24.5</td><td>4.8</td><td>3.4</td><td>0</td><td>1.3%</td></tr><tr><td>96133</td><td>1.4</td><td>1.7</td><td>1.75</td><td>24.1</td><td>4.8</td><td>3</td><td>1.0</td><td>2%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Predicted and experimental profiles for the unseeded discharge #96139. Points and errorbars refer to the experimental data, solid lines to the integrated modelling predictions. HRTS (red) and LIDAR (teal) are utilized to obtain $n_e$ and $T_e$ (figures (a) and (b)), while $T_i$ and $\Omega$ are measured by edge and core Charge Exchange (figures (c) and (d)). The predicted and interpretive regions are delimited by the vertical black lines. Since the rotation is not predicted, the solid line in the rightmost panel is the GPR fit.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Predicted and experimental profiles for the seeded discharge #96133. The notation and the layout of the plots are identical to figure 1. The rotation is not predicted and the solid line in the rightmost panel is the GPR fit. Good agreement is reached for all the predicted channels.</figcaption></figure>

QuaLiKiz a is a kinetic and electrostatic code that solves the dispersion relation to calculate the fluxes originating from microinstabilities. Specifically, it focuses on ion temperature gradient (ITG), electron temperature gradient (ETG), and trapped electron mode instabilities. QuaLiKiz utilises a quasi-linear approach, where a set of saturation rules obtained from previously performed nonlinear simulation are used on top of the linear calculations. Integrated modelling with QuaLiKiz has proven to be capable of reproducing JET profiles in various experimental conditions [47–50]. Still, extra care is needed to ensure that the cases where the QuaLiKiz assumptions break

are handled correctly. In the following, the settings chosen in presence of sawteeth, electromagnetic effect and high rotation will be described.

Modest additional transport is introduced for $\rho < 0.2$, since intermittent (1,1) magnetohydrodynamic activity (sawteeth) is not included in the simulation and the quasilinear models tend to underestimate turbulence in such region. A diffusion coefficient shaped as $H\exp(-(\rho/c_w))$ m$^2$ s$^{-1}$ is added, with $c = 0.05$ and $w = 0.15$ for all channels. $H$ is set to 0.5, 1.0 and 1.5 for $T_e$, $T_i$ and $n_e$ respectively. The same level of additional ion heat transport for $\rho < 0.2$ is set in both discharges. Since an
<table><caption>**Table 2.** Impurity composition imposed for the modelling of JET shots #96139 and #96133, normalized to the electron density. The higher tungsten content in the unseeded discharge balances out the lower temperature and lower light impurity composition, so that the radiation is similar in the two cases. The concentrations are not constant through the profiles, so average values are reported.</caption><thead><tr><th>Discharge #</th><th>Be</th><th>Ne</th><th>Ni</th><th>W</th></tr></thead><tbody><tr><td>96139</td><td>0.35%</td><td>0.05%</td><td>4.0 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>96133</td><td>1%</td><td>0.5%</td><td>4.0 × 10<sup>-4</sup></td><td>1.0 × 10<sup>-5</sup></td></tr></tbody></table>

accurate fit of the profile shapes would have been difficult, due to large errorbars and the presence of sawteeth, the consistency between the simulations was prioritized.

QuaLiKiz is an electrostatic code, but ITG driven turbulence stabilization by electromagnetic effects is important in these discharges. This effect has been seen to be non-linearly enhanced by fast ions [51]. To mimic this effect, an ad-hoc model has been proposed whereby the $R/L_{T_i}$ passed to QuaLiKiz is multiplied by the local $W_{thermal}/W_{tot}$ [50]. ETG driven mode are included in the simulations. The QuaLiKiz version applied is based on 2.8.1. However, recent work comparing QuaLiKiz and nonlinear GENE [52] highlighted an overestimation of predicted ETG fluxes in QuaLiKiz. For this reason, the electron heat flux originating from ETG is multiplied by 0.25 in all simulations compared to the QuaLiKiz 2.8.1 default.

Furthermore, this work highlighted an issue with QuaLiKiz impurity density inward convection with rotation, detailed in section 4.2, related to the breaking of the QuaLiKiz low Mach number approximation. Therefore, for the nominal simulations carried out in this paper, the impurity Mach numbers and gradient were adjusted to correspond to He masses, regardless of the actual impurity species. This version corresponds to QuaLiKiz Git commit 17c0f7cb, which in terms of physics content is identical to tagged version 2.8.1, apart from the impurity Mach number modification. Future work will focus on a more physics-based improvement to the QuaLiKiz treatment of large Mach numbers.

One of the assumptions made in [40] is that the light impurity profile is assumed to be similar to the electron profile. No experimental information about the beryllium and neon density peaking is therefore obtained using this approach. The impurity composition for the simulations shown in figures 1 and 2 is summarized in table 2. The experimental impurity composition is imposed in these runs, while the simulations with predictive impurities will be presented in section 3. Note that the different profiles obtained in predictive simulations translates in an impossibility of keeping both average $Z_{eff}$ and boundary impurity composition constant at the same time.

There is a large difference in toroidal rotation profiles between the two discharges. The value is considerably larger in the seeded case and the profiles start diverging already at the boundary conditions, hinting to a change in the rotation transport in both the SOL and the pedestal. The inversion of the rotation in the inner region is not well captured by the fit. However, QuaLiKiz does not include the impact of the rotation for $\rho < 0.4$, since the related inputs are linearly lowered to zero

in the $0.4 < \rho < 0.6$ region. This is due to the under-estimation of destabilizing parallel-velocity-gradient drive, which tends to increase with decreasing $\rho$ [53].

## 2.3. Modelling results

A prerequisite for this study is to reproduce both discharges separately in the integrated modelling framework. The agreement shown in figures 1 and 2 was considered to be sufficient for the analysis. In those simulations, and in the simulations described in this chapter, all impurities profiles were imposed from the experimental measurements. The $T_e$ and $T_i$ predictions are within the experimental uncertainties, while the density is slightly underestimated in the seeded discharge. Comparing the HRTS with the modelling of #96139, the density appears to be overestimated for $\rho < 0.5$ and underestimated for $\rho > 0.5$. Possible reasons for this disagreement will be presented in section 3. As reported in table 3, the neutron rate is slightly underestimated for #96133 (seeded), while it is overestimated for #96139 (unseeded).

Crucially, the improvement of the heat and particle core transport of the seeded pulse is captured by the modelling. $\chi_i$ drops by close to a factor two from the unseeded to the seeded discharge, from $\chi_i \sim 3.2 \pm 0.2$ (m$^2$ s$^{-1}$) to a value close to $\chi_i \sim 1.7 \pm 0.2$ (m$^2$ s$^{-1}$). These values are extracted from the simulations, at $\rho = 0.75$, and the associated error refers to the numerical fluctuations. This success is a strong indicator that the correct physics is being captured by the modelling and allows for a more detailed investigation on the causes of the improvement.

Specifically, the impact of each of the most important differences between the two discharges was assessed through a thought experiment. These differences are:

* Minor differences in the NBI and ICRH energy and total injected power

* The q profile

* The toroidal velocity profile

* The overall $Z_{eff}$ and the impurity composition

* The density pedestal

* The temperature pedestal

The strategy adopted in this paper is to start from the unseeded discharge and to gradually modify one item at a time from the list above, until the seeded discharge is obtained. The first and the last simulations of the series are the two simulations already shown in figures 1 and 2.

Various elements could in principle interact non-linearly, so the order might be important to understand the impact of every single change. Exploring all the possible combinations would however have required a large number of simulations, so the items were changed in the same order as listed above. The heating settings and the q profile are very similar in the two cases, so the first three items in the list are changed at the same time. The comparison between the profiles is shown in figure 3, while the neutron rates for all the simulations are reported in table 3.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Modelling from the step-by-step comparison between the unseeded and seeded discharges. The celeste, thicker solid and dashed lines show the baseline discharge #96139 (unseeded). The safety factor, heating and rotation profiles from #96133 (seeded) are imposed in the simulation described by the red lines. The safety factor and the heating have a very minor impact on the profiles and only the rotation is mentioned in the legend. The impurity profiles from the seeded discharge are imposed for the simulation shown in green, while the pedestal from the seeded discharge is substituted in the simulation shown by the blue lines. To emphasise the fact that this last simulation is the baseline simulation for the seeded discharge, the blue line is slightly thicker than the others. The figures on the first row show, from left to right, $n_e$, $T_e$ and $T_i$. The second row shows the toroidal rotation, $Z_{eff}$ and the safety factor $q$. The last row shows $\chi_e$ (solid) and $\chi_i$ (dashed), the electron (solid) and ion (dashed) heat fluxes and the thermal exchange between ion and electrons. In the rightmost figure on the second row positive sign represents power flowing from the electrons to the ions. Note that the impurity composition is the fixed variable, so the $Z_{eff}$ profiles are slightly different with different $T_e$ profiles.</figcaption></figure>
<table><caption>**Table 3.** Neutron rates from the step-by-step comparison between the seeded and unseeded discharges. For the first and last columns, ‘exp’ indicates the experimental neutron rates. The column designated with ‘Heating’ has the same input parameters as the simulation for #96139, but #96133 Heating settings. The simulation in the fifth column has the same inputs as the fourth, but with 96133 rotation profile. The same pattern is repeated with the impurity composition, the density pedestal and the temperature pedestal from #96133. T<sub>e</sub> and T<sub>i</sub> pedestals are changed at the same time. The values are normalized to 10<sup>15</sup> s<sup>-1</sup>.</caption><thead><tr><th>Discharge</th><th>Unseeded exp</th><th>Unseeded model</th><th>Heating & q profile</th><th>Toroidal rotation</th><th>Z<sub>eff</sub></th><th>Pedestal (seeded model)</th><th>Seeded exp</th></tr></thead><tbody><tr><td>Neutron rate (10<sup>15</sup> s<sup>-1</sup>)</td><td>4.0</td><td>5.2</td><td>5.4</td><td>5.6</td><td>5.8</td><td>10.5</td><td>10.0</td></tr></tbody></table>

<table><caption>**Table 4.** Simulation parameters corresponding to the integrated modelling, extracted at the end of the simulations when the profiles are completely relaxed and therefore stationary. The parameters for The seeded discharge are given for simulations with both interpretive and predictive impurities (imp, pred). The * next to the ion temperature gradient emphasizes the fact that the reported value is the real input to QuaLiKiz, after the application of the ad-hoc electromagnetic stabilizer. Γ<sub>E</sub> represents the E × B shearing rate.</caption><thead><tr><th>Discharge #</th><th>ρ</th><th>q</th><th>s</th><th>n<sub>e</sub> 10<sup>19</sup></th><th>T<sub>e</sub> (keV)</th><th>T<sub>i</sub> (keV)</th><th>R/L<sub>T<sub>e</sub></sub></th><th>R/L<sub>T<sub>i</sub></sub></th><th>R/L<sub>n<sub>e</sub></sub></th><th>R/L<sub>n<sub>e</sub></sub></th><th>Z<sub>eff</sub></th><th>Γ<sub>E</sub></th></tr></thead><tbody><tr><td>Int unseeded</td><td>0.75</td><td>1.99</td><td>1.87</td><td>7.36</td><td>1.00</td><td>0.94</td><td>11.1</td><td>9.68*</td><td>0.70</td><td>0.81</td><td>1.42</td><td>-0.20</td></tr><tr><td>Int seeded</td><td>0.75</td><td>2.03</td><td>1.92</td><td>7.70</td><td>1.44</td><td>1.59</td><td>12.1</td><td>9.64*</td><td>1.38</td><td>1.69</td><td>1.95</td><td>-0.39</td></tr><tr><td>Pred seeded</td><td>0.75</td><td>2.05</td><td>1.85</td><td>8.07</td><td>1.36</td><td>1.47</td><td>10.1</td><td>8.67*</td><td>2.75</td><td>9.62</td><td>1.42</td><td>-0.42</td></tr></tbody></table>

* The increase in rotation stabilizes ITG turbulence. This moderately increases n<sub>e</sub> and T<sub>i</sub>, while the increase in T<sub>e</sub> is more modest. The resulting increase in T<sub>i</sub>/T<sub>e</sub> is also stabilising for ITG. The neutron rate is slightly improved.

* The change in the impurity composition further stabilizes ITG and, to a lesser degree, ETG. This results in a further increase of T<sub>e</sub> and T<sub>i</sub>. For ITG, the fraction of the drive that is stabilized is mostly the one originating from the trapped electrons. Being its impact mostly on the electron heat flux, this increases T<sub>e</sub> more than T<sub>i</sub>. The n<sub>e</sub> profile is less peaked in the simulation with larger impurity concentration. This is consistent with the inverse dependence of density peaking on collisionality (whereas collisionality is increased by Z<sub>eff</sub>) [54, 55]. Regarding the neutron rate, the rise in the temperature and the decrease in n<sub>i</sub> balance each other out, resulting in a minor increase.

* Substituting the seeded discharge density pedestal slightly increases the overall density, since the profile peaking is unchanged and the new profile shifts almost rigidly with the boundary conditions. At the same time, the reduced NBI penetration caused by the higher density pedestal decreases the amount of fast particles in the core. The two effects balance out and the neutron rate is basically unaffected. The simulation for which only the density pedestal was modified is not shown in figure 3 for brevity.

* Changing the temperature pedestals results in a similarly rigid up-shift of the ion and electron temperatures, which lowers the collisionality even further than the unseeded discharge. This results in an increase in the density peaking and in the neutron rate.

Various parameters change between the discharges and have an impact on the turbulence. From the unseeded to the seeded, the decrease in collisionality is destabilizing, while T<sub>i</sub>/T<sub>e</sub> increases by ~10% and has a stabilizing effect. As

suggested by the integrated modelling the change in s/q is minimal, while impurities and rotation are stabilizing. Additional stabilizing mechanism for the seeded case originates from the larger normalized density gradient R/L<sub>n<sub>e</sub></sub> = R/n<sub>e</sub> dT<sub>e</sub>/dr. Lastly, due to the higher temperatures, a lower increase in the normalized temperature gradients above the threshold is enough to drive a similar heat flux in the seeded case. This is basically equivalent to a destabilization.

It is exactly because of the existence of the threshold, more than because of an exact balance of the stabilizing and destabilizing effects, that the normalized electron and ITGs are in the end similar between the unseeded and the seeded case. The parameters extracted from the integrated modelling at ρ = 0.7 are shown in table 4. Since the integrated modelling is in equilibrium and the heating is similar, those parameters do result in similar heat fluxes. Indeed, the normalized gradients are almost identical.

This stiffness also demands care in interpreting the changes in χ<sub>i</sub> and χ<sub>e</sub> shown in figure 3(g). Since a large increase in the heat flux can be caused by a small change in the gradient, for similar parameters higher fluxes correlate with higher χ<sub>s</sub>. Note in fact that, considering the normalized temperature gradient as R/L<sub>T<sub>i</sub></sub> = R/T<sub>i</sub> dT<sub>i</sub>/dr, χ<sub>i</sub> can be expressed as

$$ \chi_i = \frac{q_i R}{n T R / L_{T_i}}. \qquad (1) $$

Radiation and heating are similar in the two experiments, so the total heat flux is similar, but the heat exchange between ion and electrons varies considerably between the simulations. This effect needs to be taken into account before making conclusions about the relative importance of the stabilizing mechanisms.

To summarize, the improvement in performance, including the reduction in χ<sub>i</sub>, is mostly due stiffness and to the change in temperature pedestal height. The stabilization arising from
the different rotation and impurity composition contribute but to a lesser extent. Regarding core impurity accumulation, the simulations show that the balance between turbulence stabilization and fuel dilution leads to a rather high impurity tolerance before a loss of performance.

### 3. Sensitivities

ITG stabilization is connected to both dilution and impurity density gradient, whereby stabilization is obtained for increasing $R/L_{nimp}$ [23]. While the dilution is mostly fixed by the experimental impurity composition fixed at the boundary conditions, the peaking is more difficult to measure.

The ideal approach would be to predict the impurity transport, especially since ions are evolved by JETTO and fixing one ion might result in inconsistencies for high impurity density cases. In predictive simulation the peaking is set by the transport models. For light impurities, including neon, the neoclassical transport is sub-dominant and the impurity transport is dominated by turbulence, calculated by QuaLiKiz in our integrated modelling.

As visible in table 4, the predicted $R/L_{nNe}$ is rather large. The effect on the simulation is mostly a sharp increase in $R/L_{nNe}$, which might explain the overestimation of the density peaking in the predictive simulation. To investigate the origin of this discrepancy, it is useful to consider how the rotation is implemented in QuaLiKiz. In presence of rotation, the distribution function in QuaLiKiz can be expressed as

$$f_0 = n_0 \left( \frac{m}{2\pi T} \right)^{3/2} e^{-\frac{m}{2T}(v-U)^2}, \quad (2)$$

where $U$ is the bulk velocity of the plasma. To enable calculations, a rotation parallel to the magnetic field is assumed, $\mathbf{U} = U_{\parallel}\hat{\mathbf{b}}$ and the distribution function is expanded in the limit of small Mach number:

$$\begin{aligned} & e^{-\epsilon + \frac{m}{2}v_{\parallel}U_{\parallel} - U_{\parallel}^2} \\ & = e^{-\epsilon} \left( 1 + \frac{m}{T}v_{\parallel}U_{\parallel} + \frac{m}{2T}U_{\parallel}^2 \left( \frac{m}{T}v_{\parallel}^2 - 1 \right) \right). \end{aligned} \quad (3)$$

The linear and quadratic terms in $U_{\parallel}$ are multiplied by the mass of the ion and can become large for heavier impurities in presence of strong rotation. To avoid this problem, in nominal QuaLiKiz the mass is clamped to the mass of Beryllium for all impurities. The impurity gradient varies considerably when this parameter is modified and in this case is reduced when the mass is set to zero and increased when the true mass of the heavier impurities is used. This large impact on the impurity peaking provides an approach to perform sensitivities.

This section will study in further detail the effect of the impurities on the transport, by scanning $Z_{eff}$ and artificially manipulating the predicted density peaking. Specifically, \#96133 (seeded case) is used as the reference case and four different approaches are explored:

* Nominal JETTO-QuaLiKiz prediction

* Impurity profile imposed from the experiment as done in section 2.3, with $R/L_{nBe,Ne} = R/L_{nNe}$.

* Reduced $R/L_{nBe,Ne}$. This was carried out by entirely removing the impact of parallel flow, parallel flow shear and perpendicular flow shear on the impurity term in the dispersion relation and quasilinear flux calculations.

* Increased $R/L_{nBe,Ne}$. This was carried out by using QuaLiKiz 2.8.1, with no correction for the impurity Mach and Mach' values.

The results are summarized in table 5 and the $Z_{eff}$ scan corresponding to the decreased $R/L_{nBe,Ne}$ is shown in figure 4. $T_e$ and $T_i$ show a very consistent behaviour, increasing slightly with increasing $Z_{eff}$ in all cases. For the nominal impurity peaking most of the difference following changing $Z_{eff}$ is on the temperature channel, with ITG stabilization having little impact on the electron density profile. Even though the peaking seems closer to the experiment in the $\rho < 0.5$ region, JETTO-QuaLiKiz overestimates the total density. The large $R/L_{nNe}$ predicted in the nominal case is also not supported by the experimental data, which are not conclusive but seem to support a more reduced peaking [56]. Possible reasons for this overestimation will be explored in the next section.

When the gradient of the impurities is fixed, the change in the composition affects both the density and the temperatures. It has to be considered that the effect on the density is partially due to the fact that the simulations are not completely self-consistent. An inward pinch is in fact predicted for neon at $R/L_{Ne} = R/L_{e,0}$. This in turn translates to a larger outward main ion flux for the same electron particle flux.

Since ions are evolved in JETTO, the lower $R/L_{nNe}$ needed to drive the same flux at higher impurity concentration is partially responsible for the decrease in peaking. This effect is visible in figure 5, where the impact on the density of different impurity levels from predictive and interpretive impurity profile is shown.

The largest temperatures are reached for the simulation where $R/L_{nBe,Ne}$ are artificially increased, but the increase is not large and the neutron rate does not increase with respect to the nominal scan. Interestingly, as illustrated in figure 6, $R/L_{T_e}$ correlates with the $Z_{eff}$ absolute value more that with $R/L_{nBe,Ne}$.

### 4. Gyrokinetic analysis

#### 4.1. QuaLiKiz

In the previous sections the impurity composition was observed to have an impact on the predicted transport. A clearer picture of the mechanisms underpinning the trends observed in the integrated modelling sensitivity tests, can be drawn by analysing QuaLiKiz standalone runs. The parameters, including the impurity composition, were taken from the final stationary state of the seeded discharge \#96133 integrated modelling run with interpretive impurities, at $\rho = 0.75$. A scan of the neon density and density gradient was performed,
<table><caption>**Table 5.** Details of the Z<sub>eff</sub> scans with various assumptions on R/L<sub>nNe</sub>. Temperatures and density increase with Z<sub>eff</sub> in all cases, while the increase in the neutron rate is more contained due to the dilution. The neutron rate is largest when impurities profiles are imposed, while temperatures and densities correlate with the neon peaking in self-consistent simulations.</caption><thead><tr><th>Simulation</th><th>Z<sub>eff</sub><br>(bound)</th><th>Z<sub>eff</sub> (ave)</th><th>n<sub>e</sub> (axis),<br>10<sup>19</sup></th><th>T<sub>e</sub> (axis),<br>(KeV)</th><th>T<sub>i</sub> (axis),<br>(KeV)</th><th>Neutron<br>rate, 10<sup>15</sup></th></tr></thead><tbody><tr><td rowspan="3">Interpretive<br>R/L<sub>nNe</sub></td><td>1.4</td><td>1.55</td><td>0.94</td><td>4.7</td><td>4.6</td><td>1.04</td></tr><tr><td>1.6</td><td>1.85</td><td>0.97</td><td>5.0</td><td>4.8</td><td>1.04</td></tr><tr><td>1.8</td><td>2.15</td><td>1.01</td><td>5.2</td><td>5.0</td><td>1.04</td></tr><tr><td rowspan="3">Nominal<br>R/L<sub>nNe</sub></td><td>1.4</td><td>1.55</td><td>1.11</td><td>4.4</td><td>4.6</td><td>0.99</td></tr><tr><td>1.6</td><td>1.8</td><td>1.12</td><td>4.5</td><td>4.8</td><td>1.0</td></tr><tr><td>1.8</td><td>2.1</td><td>1.12</td><td>4.6</td><td>5.0</td><td>1.02</td></tr><tr><td rowspan="3">Decreased<br>R/L<sub>nNe</sub></td><td>1.1</td><td>1.1</td><td>1.01</td><td>4.5</td><td>4.6</td><td>0.90</td></tr><tr><td>1.4</td><td>1.4</td><td>1.02</td><td>5.0</td><td>4.9</td><td>0.94</td></tr><tr><td>1.7</td><td>1.7</td><td>1.03</td><td>5.4</td><td>5.1</td><td>0.98</td></tr><tr><td rowspan="3">Increased<br>R/L<sub>nNe</sub></td><td>1.4</td><td>1.65</td><td>1.11</td><td>4.45</td><td>4.8</td><td>0.98</td></tr><tr><td>1.6</td><td>1.95</td><td>1.12</td><td>4.7</td><td>5.0</td><td>1.0</td></tr><tr><td>1.8</td><td>2.25</td><td>1.13</td><td>4.85</td><td>5.2</td><td>1.02</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Sensitivity of the profiles to a Z<sub>eff</sub> scan for the simulations with artificially reduced R/L<sub>Ne</sub>. Electron and ion densities are shown in (a) as solid and dashed lines respectively. Electron temperature, ion temperature and Z<sub>eff</sub> profiles are shown in (b)–(d).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Electron and main ion densities for two Z<sub>eff</sub> scans. (a) shows the density profiles for the simulations where impurity transport is calculated self-consistently, while in (b) the impurity profiles are imposed. In figure (a), the Z<sub>eff</sub> indicated in the legend refers to the value imposed at the boundary. Electron and ion densities are represented by solid and dashed lines respectively. Figure (c) shows the Neon density profiles when R/L<sub>nNe</sub> = R/L<sub>nNe</sub> is imposed (dashed) and when the Neon density is predicted (solid).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Electron density, electron temperature, ion temperature and $Z_{\text{eff}}$ profiles comparing a simulation with the neon peaking predicted by nominal QuaLiKiz, in red, and a simulation with artificially reduced neon peaking in green. Average $Z_{\text{eff}}$ is similar in the two simulations, so reduced peaking also corresponds to higher $Z_{\text{eff}}$ at larger $\rho$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Fluxes as predicted by QuaLiKiz for a range of values for $Z_{\text{eff}}$ and $R/L_{n_{\text{Ne}}}$. (a)–(c) show the electron heat flux, the ion heat flux and the electron particle flux respectively. ETG is unstable for low $Z_{\text{eff}}$, while ITG is stabilized by concomitant large $Z_{\text{eff}}$ and large $R/L_{n_{\text{Ne}}}$. The fluxes are normalized in GB units.</figcaption></figure>

keeping all the other parameters fixed. The nominal point of the scan has $Z_{\text{eff}} = 1.9$ and $R/L_{n_{\text{Ne}}} = 1.7$. Rotation and the four impurity species are included.

The fluxes are expressed in Gyrobohm (GB) units, using the QuaLiKiz normalization. Heat fluxes $Q_{e,i}$ and electron particle flux $\Gamma_e$ are normalized as $Q_{e,i}_{\text{GB}} = \frac{aQ_{e,i\text{SI}}}{n_{e,i}T_{e,i}\chi_{\text{GB}}}$ and $\Gamma_{e\text{GB}} = \frac{a\Gamma_{e\text{SI}}}{n_{e,i}\chi_{\text{GB}}}$, where $a$ is the minor radius, $n_{e,i}$ the electron and ion densities and $T_{e,i}$ the electron and ion temperatures. $\chi_{\text{GB}}$ is defined as $\chi_{\text{GB}} = \frac{\sqrt{A_{i,0}m_p T_e^{1.5}}}{q_e^2 B_0^2 a}$, where $A_{i,0}$ is the mass number of the first ion, deuterium in this case, $m_p$ is the proton mass, $q_e$ is the electron charge and $B_0$ is the magnetic field.

The results are shown in figure 7. As expected, increased dilution, collisionality and $R/L_{n_{\text{Ne}}}$ are stabilizing for ITG turbulence. In the scan, the increase in neon concentration is labelled by the corresponding $Z_{\text{eff}}$, while the collisionality is changed self-consistently. The main ion density gradient goes from a maximum 1.7 to a minimum -3.4. As the impurity density gradient becomes relatively large and negative, $R/L_{n_D}$

is forced by ambipolarity to be large and positive and increasing $Z_{\text{eff}}$ is destabilizing.

At low $Z_{\text{eff}}$, a residual ETG drive carries part of the electron heat flux. This fraction of the flux is highly sensitive to $Z_{\text{eff}}$, but largely independent to the neon density gradient. It is interesting to note that at experimental levels there is still a non-negligible proportion of the electron heat flux driven by ETG. However, sensitivities in integrated modelling showed that eliminating the ETG drive does not significantly modify the profiles.

For larger impurity density gradient, increasing $Z_{\text{eff}}$ decreases the particle flux. Note that, in integrated modelling, when $R/L_{n_{\text{Ne}}}$ is larger the system increases the electron density gradient in order to preserve the same drive. In the end these two effect mostly balance out and the result is a rather constant density profile, as in figure 5(b). At large $Z_{\text{eff}}$ and low density gradient the main ion particle flux increases more quickly than the heat flux, leading to a flatter density profile, in accord with figure 4(a).
<table><caption>**Table 6.** Parameters used for the comparison between GENE and QuaLiKiz. *x* is the local normalized radius, *x* = *a*/R<sub>min</sub>, with the minor radius R<sub>min</sub> = 0.89 (m).</caption><thead><tr><th>x</th><th>ρ</th><th>q</th><th>s</th><th>n<sub>e</sub> 10<sup>19</sup></th><th>T<sub>e</sub> (keV)</th><th>T<sub>i</sub> (keV)</th><th>R/L<sub>T<sub>e</sub></sub></th><th>R/L<sub>T<sub>i</sub></sub></th><th>R/L<sub>n<sub>e</sub></sub></th><th>R/L<sub>η<sub>Ne</sub></sub></th><th>Z<sub>eff</sub></th><th>Γ<sub>E</sub></th></tr></thead><tbody><tr><td>0.76</td><td>0.7</td><td>1.74</td><td>1.69</td><td>8.27</td><td>1.62</td><td>1.66</td><td>10.53</td><td>9.84</td><td>2.01</td><td>1.62</td><td>2.15</td><td>-0.41</td></tr></tbody></table>

<table><caption>**Table 7.** Grid parameters used for the GENE simulation. *n*<sub>x</sub>, *n*<sub>y</sub>, *n*<sub>z</sub>, *n*<sub>v</sub> and *n*<sub>w</sub> represent the number of grid points respectively for the radial, bi-normal, parallel, *v*<sub>||</sub> and magnetic moment dimensions. The large *n*<sub>x</sub> is necessary to resolve the Floquet cycles. *L*<sub>x</sub>, *L*<sub>v</sub> and *L*<sub>w</sub> are the extension of the simulation box in the radial, *v*<sub>||</sub> and magnetic moment directions, normalized to the inverse gyroradius, thermal velocity and *T*<sub>e</sub>/B respectively.</caption><thead><tr><th>n<sub>x</sub></th><th>n<sub>z</sub></th><th>n<sub>v</sub></th><th>n<sub>w</sub></th><th>L<sub>x</sub></th><th>L<sub>v</sub></th><th>L<sub>w</sub></th></tr></thead><tbody><tr><td>151</td><td>32</td><td>64</td><td>24</td><td>40.0</td><td>3.0</td><td>9.0</td></tr></tbody></table>

<table><caption>**Table 8.** Comparison between GENE and QuaLiKiz growthrates and frequencies, expressed in GENE gyroBohm units. Agreement in γ<sub>GB</sub> is obtained in the last row by decreasing the QuaLiKiz input R/L<sub>T<sub>i</sub></sub> by 20%.</caption><thead><tr><th></th><th>γ<sub>GB</sub></th><th>ω<sub>GB</sub></th></tr></thead><tbody><tr><td>QuaLiKiz nominal</td><td>0.574</td><td>1.27</td></tr><tr><td>QuaLiKiz -20%R/L<sub>T<sub>i</sub></sub></td><td>0.306</td><td>1.19</td></tr><tr><td>GENE</td><td>0.323</td><td>0.926</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Transport coefficients for QuaLiKiz and GENE with rotation and ExB shear set to zero. Growthrates and v<sub>Ne</sub>/D<sub>Ne</sub>, which gives information about the peaking, are reported in (a) and (b). The diffusion D<sub>Ne</sub> and the pinch coefficient v<sub>Ne</sub> are reported in (c) and (d), divided by the main ion heat flux. k<sub>θ</sub>ρ<sub>s</sub> is scanned up to the last unstable growthrate at k<sub>θ</sub>ρ<sub>s</sub>. The solid blue and green lines are QuaLiKiz and GENE with identical parameters, while for the simulations represented by the blue dashed line the QuaLiKiz input R/L<sub>T<sub>i</sub></sub> is reduced by 15%. All values are reported with the GyroBohm normalization.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9. Transport coefficients for QuaLiKiz and GENE with $\Gamma_E$ set to zero and without $E \times B$ shear respectively. The structure of the plots is the same as in figure 8. The solid blue and green lines are QuaLiKiz and GENE with identical parameters, while for the simulations represented by the blue dashed line the QuaLiKiz input $R/L_T$ is increased by 5%. All values are reported with the GyroBohm normalization. With respect to the GENE prediction, $D_{Ne}/\chi$ is underestimated by QuaLiKiz.</figcaption></figure>

## 4.2. GENE

From the analysis presented above, the importance of QuaLiKiz being able to correctly predict the impurity transport is evident. This section will focus on the effect of the rotation on the impurity profiles, since as shown in figure 6, section 3, it was found to be a strong actuator. Other phenomena not captured in QuaLiKiz can impact the impurity transport, as for example finite $\beta$ effects [57] and fast ions [58], but since they are not expected to be relevant in this case such effects will not be included in this work.

The QuaLiKiz predictions are compared with the output from the higher fidelity gyrokinetic code GENE. Only the linear physics is validated and direct comparison to nonlinear simulations, is left for future work. For simplicity, only Neon impurities are considered. $\beta$ is set to zero and the $s-\alpha$ geometry model is used. The collisions are treated by a linearized Landau-Boltzmann operator. The input parameters are extracted from integrated modelling and are summarized in table 6. The grid parameters for the simulations are listed in table 7. Henceforth, the growthrates will be presented using the GENE GB normalization employed in the simulation.

GENE and QuaLiKiz have already been extensively compared [59–61], including comparisons of the stabilization

from increasing $Z_{eff}$ [62]. A basic level of agreement in this corner of the parameter space is confirmed here by comparing the stabilization that originates from the addition of Neon. For $k_\theta \rho_s = 0.2$, which is close to the most unstable mode, $\gamma_{GB}$ drops by 25% in GENE and by 20% in QuaLiKiz. Here $k_\theta$ is the poloidal wavenumber and $\rho_s = \frac{\sqrt{T_e m_i}}{q_e B}$ the main ion gyroradius. A reduction in $R/L_T$ of 15% is sufficient to match the growthrates. The results are summarized in table 8.

$D_s$ and $V_s$ are calculated by QuaLiKiz, but are not provided directly by linear GENE. While the absolute values of the saturated fluxes cannot be obtained without the nonlinear terms, it is still possible to extract $D_s/\chi$ and $V_s/\chi$, since these transport coefficients depend only on the ratio of the fluxes. To extract $D$ and $V$ in the GENE linear simulations, an extra trace neon species was added with zero gradient, and the neon transport coefficients were calculated using

$$V_{Ne} = \Gamma_{\text{trace}}/n_{\text{trace}}, \quad (4)$$

$$D_{Ne} = (\Gamma_{\text{main}} - n_{\text{main}}V_{Ne})(R_0/(n_{\text{main}}R/L_{n,Ne})). \quad (5)$$

$D_{Ne}/\chi$ and $V_{Ne}/\chi$ are calculated for $k_\theta \rho_s = 0.1, 0.2, 0.3, 0.4$ and the results are compared with the
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Transport coefficients for QuaLiKiz and GENE with the $E \times B$ shear divided by 4 with respect to the nominal value. The structure of the plots is the same as in figure 8. The green, blue and red lines represent GENE, QuaLiKiz and QuaLiKiz with the effect of the rotation on the neon term in the dispersion relation and quasilinear flux calculations removed. The input $R/L_{T_i}$ for the two QuaLiKiz versions is reduced by 4% and 7% respectively in order to match the GENE growthrates. All values are reported with the GyroBohm normalization. With respect to the GENE prediction, $D_{Ne}/\chi$ and $v_{Ne}/D_{Ne}$ are underestimated by QuaLiKiz. The inward peaking is overestimated for nominal QuaLiKiz and underestimated for the customized version.</figcaption></figure>

QuaLiKiz predictions. Agreement in the growthrates is reached by decreasing the $R/L_{T_i}$ in QuaLiKiz by 20%. Considering the two cases with the same drive, both $D_{Ne}/\chi$ and $v_{Ne}/D_{Ne}$ are slightly overestimated and the $v_{Ne}/D_{Ne}$ ratio is very similar. Growthrates and transport coefficients for the rotation-less case are reported in figure 8.

The effect of rotation on the impurity peaking is then analysed. First, the comparison is performed without $E \times B$ shear in GENE and with $\gamma_E$ set to zero in QuaLiKiz. This simplifies the analysis of the GENE simulations and the results are reported in figure 9.

QuaLiKiz and Gene are in good agreement with regards to the particle pinch, but QuaLiKiz underestimates the diffusion coefficient by a factor ~2. This is not due to the difference in growthrates, since the same drive was retrieved by changing the input ITG of just 5%. Interestingly, the $v/D$ predicted by Gene does not change significantly from the rotation-less case. The disagreement is not related to the impurity mass in the low Mach number approximation and it possibly hints to an overestimation of the impact of the parallel velocity gradient in QuaLiKiz. The identification of the exact reason for the

low $D_{Ne}$ is left for future work. This difference explains the overestimation of the impurity peaking in integrated modeling with nominal QuaLiKiz.

Lastly, the effect of the $E \times B$ shear on the particle transport coefficients is studied. The GENE simulations including this effect were analyzed with the methodology described in appendix. Unfortunately, such methodology proved not to be suitable at nominal rotation. Within the scope of code comparison between GENE and QuaLiKiz, the predicted transport coefficients were compared at lower $E \times B$ shear. The $E \times B$ shear was therefore reduced by a factor 4 with respect to the experimental value. The differences in the prediction of growthrates and transport coefficient ratios between GENE and QuaLiKiz at reduced shear can be observed comparing the green and blue lines in figure 10. A 4% reduction of the input $R/L_{T_i}$ passed to QuaLiKiz was sufficient to obtain similar growthrates. Since this had little to no effect on the neon transport coefficients ratios, the simulations with the same drive are shown for brevity.

Some differences between the two models are evident. At larger $k_\theta \rho_s$, QuaLiKiz is stabilised faster than GENE. The
pinch coefficient corresponding to the largest growth rates is underestimated by QuaLiKiz by only a few percent, but there is a factor ~2 difference in the diffusivity. In integrated modelling, this corresponds to an overestimation of the peaking of the impurity profile by QuaLiKiz.

The impact of rotation on the impurity transport in QuaLiKiz with the effect of rotation on the impurities removed in the dispersion relation and quasilinear flux calculations is shown by the solid red line. In this case, a 7% reduction of the input $R/L_{T_i}$ was sufficient to match the growthrates and again only the simulations with the same drive are shown. $D_{Ne}$ is similar, but the pinch coefficient absolute value is greatly reduced and the corresponding peaking in the integrated modelling is in this case underestimated.

The results at lower rotation point to a discrepancy between GENE and QuaLiKiz with respect to the impact of rotation on impurity transport, consistent with the exaggerated Neon peaking predicted by QuaLiKiz with nominal rotation.

## 5. Conclusions

Integrated modelling of plasma discharges including multiple impurities was performed. Specifically, the role of neon in modifying the core turbulent transport was analysed. This was motivated by the experimental observation of an increase in confinement with neon seeding, for both pedestal and core. A model to understand the changes in the pedestal, where most of the improvement in confinement and neutron rate originates, is not available. Future work might focus on possible mechanisms impacting the pedestal performance, in particular:

* In the unseeded case the maximum pressure gradient was closer to the separatrix, but the pedestal was narrower. The two effects are respectively destabilizing and stabilizing.

* A decrease of $T_c$ at the separatrix, due to the radiation and increased penetration of recycling neutrals, is associated with an inward shift of the pressure profile, which is a potential stabilization mechanism [37]. The experimental assessment of its importance is however difficult, since the position of the separatrix is uncertain.

* The increase in $Z_{eff}$ and dilution might also stabilize pedestal turbulence, which would increase the pedestal height. Pedestal turbulence stabilization due to higher $Z_{eff}$ is indeed expected according to gyrokinetic calculations [63] and simulations [64–66].

Regarding the core, the transport improvement is identified to be due to an increase in rotation and from the impurity-induced stabilization of ITG modes. The higher ion temperature balances the dilution, leaving the neutron rate mostly unchanged. $T_c$ and $T_i$ are well predicted, while $n_e$ is overestimated when the impurity transport is predicted and underestimated when impurities are imposed. This is linked to an overestimation of the light impurity density peaking by QuaLiKiz.

The toroidal rotation was identified as a strong driver of the neon inward flux. The QuaLiKiz prediction was compared to higher fidelity GENE calculations, showing a lower increase in impurity peaking with rotation in the higher fidelity code. Future QuaLiKiz development will address this difference. Comparison with other quasilinear models, such as TGLF [67], could at the same time aid with the interpretation of the experiment and provide an alternative prediction that can be used to identify broken assumptions and possible improvements for both codes.

While more work needs to be done to understand the behaviour of the pedestal and the extrapolability of the observed improvement, the results of this paper is optimistic regarding the maximum tolerable light impurity content that can be injected in reactors to cool the divertor region.

## Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. A previous version of this paper appeared in the thesis of the first author M Marin, published in September 2021 with the title *Integrated modelling of multiple-ion tokamak discharges*.

## Appendix

In order to extract the impurity transport coefficients from linear GENE simulations including the $E \times B$ shear, the following needs to be considered.

The $E \times B$ shear in GENE is implemented as a shift in $k_x$. Every time-step $k_x$ is shifted by $\Delta k_x = \frac{\Gamma \gamma_E}{R_0 c_{ref}}$, with $t$ being the time-step in GENE, $R_0$ the major radius and $c_{ref} = \sqrt{\frac{q_e T_{ref}}{m_{ref}}}$. Since the grid is discrete, the shift is implemented only when it accumulates for a value larger than $\Delta k_x \equiv 2\pi/L_x$.

The $k_x$ spectra for $k_y = 0.3$ as predicted by GENE is shown in figure 11. The maximum in the instability is found at $\Delta k_x = \pm 0.35$ and decreases sharply for higher shifts. While the mode traverses the whole spectra in the simulations including the $E \times B$ shear, the growthrate rises and falls. This behaviour creates cyclical oscillations called Floquet modes. The impurity transport coefficients vary during the cycles, with both $D_{Ne}$ and $V_{Ne}$ increasing with $\Delta k_x$. $D/V$ is non monotonic and has a maximum around 0.5.

The rotation is therefore included utilizing the $\tau_{AC}$ approach, as presented in [68]. This method originates from the observation that the relevant timescale to average physical quantities is the nonlinear decorrelation time. If the mode does not maintain coherence for long enough to traverse the full Floquet cycle, it is reasonable to expect that nonlinearly
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Details of the $k_x$ spectra for $k_y = 0.3$ as predicted by GENE. The structure of the plots is the same as in figure 8. Half of the period is shown.</figcaption></figure>

higher $k_x$ would be less important. The decorrelation time can be estimated by $1/\gamma_{k_y}$ and is usually minimum around $k_x = 0$. Taking these considerations into account, the transport coefficients are calculated by averaging their values over the nonlinear decorrelation time.

The following procedure is therefore repeated for the four different binormal wavenumbers, $k_\theta \rho_s = 0.1, 0.2, 0.3, 0.4$:

* Linear GENE is run with the appropriate $E \times B$ shear until a few Floquet cycles are computed.

* The growthrates for each time-step are calculated as $\gamma = \ln(\frac{n_{ij}(t+\Delta t)}{n_{ij}(t)})/\Delta t$.

* The discrete steps in $k_x$ result in large growthrates when the distribution function is shifted. After a few time-steps, a new configuration is reached and the growthrate returns to the pre-shift levels. To allow the identification of the maxima of the Floquet cycles, a smoothing over the discrete shift timescale is performed. The new growthrates are calculated as $\gamma = \ln(\frac{n_{ij}(t+n\Delta t)}{n_{ij}(t)})/(n\Delta t)$. n here is defined as $n = T/\Delta t$, where T is the time between two consequent shifts.

* The maxima in the growth rates are identified, now corresponding to the most unstable $k_x$, usually close to $k_x \sim 0$. Then, for each cycle:

* A time window satisfying the relation $|1 - \gamma_{k_y} \Delta \tau| < \delta$ is identified. Here $\delta$ is a small parameter and $\tau$ is the width of the time window where $\gamma_{k_y}$ is averaged. The centre of the window is kept fixed at the time corresponding to the maximum growthrate and the width is increased symmetrically until the relation is satisfied.

* Inside of this window, for each time-step, $V_{Ne}/\chi$ and $D_{Ne}/\chi$ are calculated.

* To obtain the final $V_{Ne}/\chi$ and $D_{Ne}/\chi$, the entirety of the corresponding values inside the identified intervals are averaged.

This last step is justified when neither coefficient changes significantly inside the intervals under consideration.

The results for $k_y = 0.2$ are shown in figure 12. The growthrates are well converged for all the cycles and $V_{Ne}$ and $D_{Ne}$ are fairly constant over the intervals considered.

As shown in figure 13, with $k_\theta \rho_s = 0.2$, the particle transport coefficients for nominal rotation are not constant over the averaged time windows and vary considerably during the Floquet cycle. Part of the reason originates from the $E \times B$ stabilization itself, which results in lower $\gamma_k$. This in turn translates into a larger $\tau$, thus widening the range where the values are averaged. Nonlinear simulations might be needed for comparison at large rotation values and will be left for future work.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Visualization of the calculations performed to obtain the neon transport coefficients for $k_{\theta} \rho_{s} = 0.2$. Plots for only one poloidal wavenumber at low $E \times B$ shear simulations are shown, the other are similar and are not shown for brevity. Plot (a) reports the growthrates extracted during multiple Floquet cycles in blue and the fit performed to locate the maxima in green. The growthrates here have already been smoothed according to item 3 of the GENE methodology list. The interval centered around the maxima and satisfying $|1 - \gamma_{k} \Delta \tau| < \delta$ are indicated by the black vertical lines. $V_{Ne}/Q_{i}$ and $D_{Ne}/Q_{i}$ are shown for all the time steps of the simulation in (b) and only inside the selected intervals in (c).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Visualization of the calculations performed to obtain the neon transport coefficients for $k_{\theta} \rho_{s} = 0.2$ at nominal $E \times B$ shear. The structure of the plots is the same as in figure 12. The black intervals satisfying $|1 - \gamma_{k} \Delta \tau| < \delta$ are visibly larger, resulting in non constant transport coefficients in (c).</figcaption></figure>

## ORCID iDs

J. Citrin https://orcid.org/0000-0001-8007-5501

C. Bourdelle https://orcid.org/0000-0002-4096-8978

L. Garzotti https://orcid.org/0000-0002-3796-9814

A. Ho https://orcid.org/0000-0001-5107-3531

## References

[1] Pitts R.A. et al 2019 Physics basis for the first ITER tungsten divertor *Nucl. Mater. Energy* **20** 100696

[2] Eich T. et al 2013 Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER *Nucl. Fusion* **53** 093031

[3] Asakura N. et al 2017 Studies of power exhaust and divertor design for a 1.5 GW-level fusion power DEMO *Nucl. Fusion* **57** 126050

[4] Wenninger R.P. et al 2014 DEMO divertor limitations during and in between ELMs *Nucl. Fusion* **54** 114003

[5] Kallenbach A. et al 2021 Developments towards an ELM-free pedestal radiative cooling scenario using noble gas seeding in ASDEX Upgrade *Nucl. Fusion* **61** 016002

[6] Giroud C. et al 2012 Integration of a radiative divertor for heat load control into JET high triangularity ELMy H-mode plasmas *Nucl. Fusion* **52** 063022

[7] Jaervinen A.E. et al 2016 Impact of divertor geometry on radiative divertor performance in JET H-mode plasmas *Plasma Phys. Control. Fusion* **58** 045011

[8] Kubo H. et al 2001 High radiation and high density experiments in JT-60U *Nucl. Fusion* **41** 227-33

[9] Jackson G.L. et al 2002 Effects of impurity seeding in DIII-D radiating mantle discharges *Nucl. Fusion* **42** 28-41

[10] Tardini G. et al 2013 Core transport analysis of nitrogen seeded H-mode discharges in the ASDEX Upgrade *Plasma Phys. Control. Fusion* **55** 015010

[11] Giroud C. et al 2013 Impact of nitrogen seeding on confinement and power load control of a high-triangularity JET ELMy H-mode plasma with a metal wall *Nucl. Fusion* **53** 113025

[12] Dunne M.G. et al 2017 Global performance enhancements via pedestal optimisation on ASDEX Upgrade *Plasma Phys. Control. Fusion* **59** 025010

[13] Casali L., Fable E., Dux R. and Ryter F. 2018 Modelling of nitrogen seeding experiments in the ASDEX Upgrade tokamak *Phys. Plasmas* **25** 032506

[14] Pütterich T. et al 2013 Observations on the W-transport in the core plasma of JET and ASDEX Upgrade *Plasma Phys. Control. Fusion* **55** 124036

[15] Drenik A. et al 2019 Evolution of nitrogen concentration and ammonia production in N2-seeded H-mode discharges at ASDEX Upgrade *Nucl. Fusion* **59** 046010

[16] Telesca G., Ivanova-Stanik I., Brezinsek S., Czarnecka A., Drewelow P., Giroud C., Huber A., Wiesen S., Wischmeier M. and Zagórski R. 2017 Simulation of JET ITER-Like Wall pulses at high neon seeding rate *Nucl. Fusion* **57** 126021

[17] Beurskens M.N. et al 2008 Pedestal and ELM response to impurity seeding in JET advanced scenario plasmas *Nucl. Fusion* **48** 095004

[18] Ongena J. et al 1996 Confinement transitions with radiation cooling in TEXTOR-94 *Plasma Phys. Control. Fusion* **38** 279-88

[19] Messiaen A.M. et al 1996 High confinement and high density with stationary plasma energy and strong edge radiation in the TEXTOR-94 Tokamak *Phys. Rev. Lett.* **77** 2487-90

[20] McKee G.R. et al 2000 Impurity-induced turbulence suppression and reduced transport in the DIII-D tokamak *Phys. Plasmas* **7** 1870-7

[21] Glöggler S. et al 2019 Characterisation of highly radiating neon seeded plasmas in JET-ILW *Nucl. Fusion* **59** 126031

[22] Angioni C. 2022 Impurity transport in tokamak plasmas, theory, modelling and comparison with experiments *Nucl. Fusion* **63** 073001

[23] Bonanomi N., Mantica P., Citrin J., Giroud C., Lerche E., Sozzi C., Taylor D., Tsalas M. and Van Eester D. 2018
Effects of nitrogen seeding on core ion thermal transport in JET ILW L-mode plasmas Nucl. Fusion **58** 026028

[24] Moradi S., Bourdelle C., Tokar M.Z., Litaudon X., Imbeaux F., Corre Y., Monier-Garbet P., Kalupin D. and Weyssow B. 2012 Modeling of energy confinement improvement in high power JET discharges with neon seeding *Plasma Phys. Control. Fusion* **54** 015004

[25] Li J., Wang Z.X., Dong J.Q., Han M.K., Shen Y., Xiao Y. and Du H.R. 2019 Impurity effects on ion temperature gradient driven multiple modes in transport barriers Nucl. Fusion **59** 076013

[26] Sciortino F. et al 2020 Inference of experimental radial impurity transport on alcator C-Mod: Bayesian parameter estimation and model selection Nucl. Fusion **60** 126014

[27] Meneghini O., Garud S., Lyons B., McClenaghan J., Imai C., Grierson B., Smith S., Staebler G.M., Snyder P. and Candy J. 2021 Neural-network accelerated coupled core-pedestal simulations with self-consistent transport of impurities and compatible with ITER IMAS Nucl. Fusion **61** 026006

[28] Koechl F. et al 2018 Optimising the ITER 15MA DT baseline scenario by exploiting self-consistent-consistent free-boundary core-edge-sol workflow in IMAS IAEA Fusion Energy Conf. (Gandhinagar, India, 22–27 October 2018) pp 7–25

[29] Luda T., Angioni C., Dunne M., Fable E.E., Kallenbach A., Bonanomi N., Schneider P.A., Siccinio M. and Tardini G. (The ASDEX Upgrade Team and The EUROfusion MST1 Team) 2020 Integrated modeling of ASDEX Upgrade plasmas combining core, pedestal and scrape-off layer physics Nucl. Fusion **60** 036023

[30] Citrin J. et al 2017 Tractable flux-driven temperature, density and rotation profile evolution with the quasilinear gyrokinetic transport model QuaLiKiz *Plasma Phys. Control. Fusion* **59** 124005

[31] Jenko F., Dorland W., Kotschenreuther M. and Rogers B.N. 2000 Electron temperature gradient driven turbulence *Phys. Plasmas* **7** 1904–10

[32] Ho A., Citrin J., Auriemma F., Bourdelle C., Casson F., Kim H.-T., Manas P., Szepesi G. and Weisen H. 2019 Application of Gaussian process regression to plasma turbulent transport model validation via integrated modelling Nucl. Fusion **59** 056007

[33] Lao L.L., John H.S., Stambaugh R.D., Kellman A.G. and Pfeiffer W. 1985 Reconstruction of current profile parameters and plasma shapes in tokamaks Nucl. Fusion **25** 1611

[34] Szepesi G., Appel L.C., Luna E.D., Gelfusa M., Hawkes N.C., Gaudio P. and Sertoli M. 2020 Advanced equilibrium reconstruction for JET with EFIT++ 47th EPS Conf. on Plasma Physics Online Conf. (21–25 June 2021) p 1037

[35] Stangeby P.C. 2000 *The Plasma Boundary of Magnetic Fusion Devices* (London: Taylor and Francis)

[36] Frassinetti L. et al 2019 Role of the pedestal position on the pedestal performance in AUG, JET-ILW and TCV and implications for ITER Nucl. Fusion **59** 076038

[37] Saarelma S., Järvinen A., Beurskens M., Challis C., Frassinetti L., Giroud C., Groth M., Leyland M., Maggi C. and Simpson J. 2015 The effects of impurities and core pressure on pedestal stability in Joint European Torus (JET) Phys. Plasmas **22** 2–7

[38] Romanelli M. et al 2014 JINTRAC: a system of codes for integrated simulation of Tokamak scenarios *Plasma Fusion Res.* **9** 3403023

[39] Romanelli G.C.M. 1988 *Jetto Manual* (Italy: JET Joint Undertaking) p 84

[40] Sertoli M., Carvalho P.J., Giroud C. and Menmuir S. 2019 Measuring the plasma composition in tokamaks with

metallic plasma-facing components *J. Plasma Phys.* **85** 22–27

[41] Houlberg W.A., Shaing K.C., Hirshman S.P. and Zarnstorff M.C. 1997 Bootstrap current and neoclassical transport in tokamaks of arbitrary collisionality and aspect ratio *Phys. Plasmas* **4** 3230

[42] Challis C.D., Cordey J.G., Hamnén H., Stubberfield P.M., Christiansen J.P., Lazzaro E., Muir D.G., Stork D. and Thompson E. 1989 Non-inductively driven currents in JET Nucl. Fusion **29** 563–70

[43] Eriksson L.G., Hellsten T. and Willen U. 1993 Comparison of time dependent simulations with experiments in ion cyclotron heated plasmas Nucl. Fusion **33** 1037–48

[44] Lauro-Taroni L., Alper B., Giannella R., Marcus F., Smeulders P., Von Hellermann M., Lawson K. and Mattioli M. 1994 Impurity transport of high performance discharges in JET 21st EPS Conf. on Controlled Fusion and Plasma Physics (Montpellier, France, 27 June–1 July 1994)

[45] Tamor S. 1981 ANTIC: A code for Calculation of Neutral Transport in Cylindrical Plasmas *J. Comput. Phys.* **40** 104119

[46] Oberparleiter M., Jenko F., Told D., Doerk H. and Görler T. 2016 Interaction between neoclassical effects and ion temperature gradient turbulence in gradient- and flux-driven gyrokinetic simulations *Phys. Plasmas* **23** 1–9

[47] Baiocchi B., Bourdelle C., Angioni C., Imbeaux F., Loarte A. and Maslov M. 2015 Transport analysis and modelling of the evolution of hollow density profiles plasmas in JET and implication for ITER Nucl. Fusion **55** 123001

[48] Breton S. et al 2018 First principle integrated modeling of multi-channel transport including Tungsten in JET Nucl. Fusion **58** 096003

[49] Linder O., Citrin J., Hogeweij G., Angioni C., Bourdelle C., Casson F., Fable E., Ho A., Koechl F. and Sertoli M. 2019 Flux-driven integrated modelling of main ion pressure and trace tungsten transport in ASDEX Upgrade Nucl. Fusion **59** 016003

[50] Casson F.J. et al (JET Contributors) 2020 Predictive multi-channel flux-driven modelling to optimise ICRH tungsten control and fusion performance in JET Nucl. Fusion **60** 066029

[51] Citrin J., Jenko F., Mantica P., Told D., Bourdelle C., Garcia J., Haverkort J.W., Hogeweij G.M., Johnson T. and Pueschel M.J. 2013 Nonlinear stabilization of tokamak microturbulence by fast ions *Phys. Rev. Lett.* **111** 155001

[52] Citrin J. et al (JET Contributors) 2021 Integrated modelling and multiscale gyrokinetic validation study of ETG turbulence in a JET hybrid H-mode scenario Nucl. Fusion **62** 086025

[53] Citrin J. et al 2017 Tractable flux-driven temperature, density and rotation profile evolution with the quasilinear gyrokinetic transport model QuaLiKiz *Plasma Phys. Control. Fusion* **59** 124005

[54] Angioni C. et al 2007 Scaling of density peaking in H-mode plasmas based on a combined database of AUG and JET observations Nucl. Fusion **47** 1326–35

[55] Angioni C., Fable E., Greenwald M., Maslov M., Peeters A.G., Takenaga H. and Weisen H. 2009 Particle transport in tokamak plasmas, theory and experiment *Plasma Phys. Control. Fusion* **51** 124017

[56] Giroud C. 2022 in preparation

[57] Hein T., Angioni C., Fable E. and Candy J. 2010 Gyrokinetic study of the role of $\beta$ on electron particle transport in tokamaks *Phys. Plasmas* **17** 102309

[58] Manas P., Kappatou A., Angioni C. and McDermott R.M. 2020 Light impurity transport in tokamaks: On the impact of neutral beam fast ions Nucl. Fusion **60** 056005

[59] Bourdelle C., Garbet X., Imbeaux F., Casati A., Dubuit N., Guirlet R. and Parisot T. 2007 A new gyrokinetic
quasilinear transport model applied to particle transport in tokamak plasmas *Phys. Plasmas* **14** 112501

[60] Casati A. et al 2009 Validating a quasi-linear transport model versus nonlinear simulations *Nucl. Fusion* **49** 085012

[61] Cottier P., Bourdelle C., Camenen Y., Gürcan O.D., Casson F.J., Garbet X., Hennequin P. and Tala T. 2014 Angular momentum transport modeling: achievements of a gyrokinetic quasi-linear approach *Plasma Phys. Control. Fusion* **56** 015011

[62] Bourdelle C., Garbet X., Hoang G., Ongena J. and Budny R. 2002 Stability analysis of improved confinement discharges: internal transport barriers in Tore Supra and radiative improved mode in TEXTOR *Nucl. Fusion* **42** 892–902

[63] Dominguez R.R. and Rosenbluth M.N. 1989 Local kinetic stability analysis of the ion temperature gradient mode *Nucl. Fusion Lett.* **29** 844

[64] Kotschenreuther M., Hatch D.R., Mahajan S., Valanju P., Zheng L. and Liu X. 2017 Pedestal transport in H-mode plasmas for fusion gain *Nucl. Fusion* **57** 064001

[65] Hatch D.R., Kotschenreuther M., Mahajan S., Valanju P. and Liu X. 2017 A gyrokinetic perspective on the JET-ILW pedestal *Nucl. Fusion* **57** 036020

[66] Hatch D.R. et al 2019 Direct gyrokinetic comparison of pedestal transport in JET with carbon and ITER-like walls *Nucl. Fusion* **59** 086056

[67] Staebler G.M., Candy J., Howard N.T. and Holland C. 2016 The role of zonal flows in the saturation of multi-scale gyrokinetic turbulence *Phys. Plasmas* **23** 062518

[68] Dagnelie V.I., Citrin J., Jenko F., Pueschel M.J., Görler T., Told D. and Doerk H. 2019 Growth rates of ITG modes in the presence of flow shear *Phys. Plasmas* **26** 046010