# Full radius integrated modelling of ohmic ramp-up at TCV including self consistent density prediction

M. Marin, Yann Camenen, C. Bourdelle, F.J. Casson, R. Coosemans, L.

Garzotti, P. Maget, P. Manas, A. Najlaoui, O. Sauter

► To cite this version:

M. Marin, Yann Camenen, C. Bourdelle, F.J. Casson, R. Coosemans, et al. Full radius integrated modelling of ohmic ramp-up at TCV including self consistent density prediction. Nuclear Fusion, 2025, 65, pp.036015. 10.1088/1741-4326/adb169. hal-05372865

HAL Id: hal-05372865

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

Submitted on 19 Nov 2025

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 Attribution 4.0 International License
PAPER • OPEN ACCESS

# Full radius integrated modelling of ohmic ramp-up at TCV including self consistent density prediction

To cite this article: M. Marin et al 2025 *Nucl. Fusion* 65 036015

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.

Mass spectrometers for vacuum, gas, plasma and surface science

HIDEN
ANALYTICAL

# Ultra-high Resolution Mass Spectrometers for the Study of Hydrogen Isotopes and Applications in Nuclear Fusion Research

## DLS Series

* **Unique** Dual Mass range / Zone H functionality

* For the measurement of overlapping species

* He/D2, CH2D2/H20, Ne/D20

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

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

## HAL 101X

* Monitoring, diagnostics and analysis applications in tokamak and torus operations

* Unique design avoids all radiation shielding requirements

* Featuring TIMS mode for real-time quantification of hydrogen and helium isotopes

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

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

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

www.HidenAnalytical.com

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

info@hiden.co.uk
# Full radius integrated modelling of ohmic ramp-up at TCV including self consistent density prediction

M. Marin<fnref n="1" />,*, Y. Camenen<fnref n="2" />, C. Bourdelle<fnref n="3" />, F.J. Casson<fnref n="4" />, R. Coosemans<fnref n="1" />, L. Garzotti<fnref n="4" />, P. Maget<fnref n="3" />, P. Manas<fnref n="2" />, A. Najlaoui<fnref n="2" />,<fnref n="5" />, O. Sauter<fnref n="1" /> and the TCV Teamª

E-mail: michele.marin@epfl.ch

Received 19 June 2024, revised 24 January 2025

Accepted for publication 3 February 2025

Published 14 February 2025

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

CrossMark

## Abstract

The ramp-up is a critical phase in the operations of a Tokamak, during which engineering and physics aspects must be taken into account to ensure stability, minimize flux consumption and avoid disruptions. Predicting ramp-up phases faces challenges such as nonlinearity, uncertainty on boundary and initial conditions and changes in the magnetic equilibrium. Our work uses the High-Fidelity Pulse Simulator (HFPS), a Python workflow based on JINTRAC. The input and output are in machine and code generic IMAS data format. The HFPS predicts the evolution of the current, temperature, main ion density and impurity density up to the separatrix. The self-consistent prediction of the density during the ramp-up represents the main element of novelty in this work. To this end, a closed feedback loop is set to match experimental line averaged density. QuaLiKiz (Citrin et al 2017 *Plasma Phys. Control. Fusion* **59** 124005), TGLF (Staebler et al 2016 *Phys. Plasmas* **23** 062518) and FRANTIC (Tamor 1981 *J. Comput. Phys.* **40** 104119) are used to calculate turbulent fluxes and the source of neutrals respectively. QuaLiKiz and TGLF predict a transition from Trapped Electron Mode early in the discharge to Ion Temperature Gradient dominated turbulence. The results are compared to higher fidelity simulations with GKW, which show qualitative agreement. Good general agreement is reached between integrated modelling and experimental data, quantified by proposed measures of agreement. A large set of sensitivities to modelling choices and initial and boundary conditions is performed on four different discharges, to assess the robustness of the approach.

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

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

<footnote n="2">² CNRS, Aix-Marseille University, PIIM UMR7345, Marseille, France</footnote>

<footnote n="3">³ CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France</footnote>

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

<footnote n="5">⁵ ITER Organization, Route de Vinon-sur Verdon, 13067 St. Paul Lez Durance, France</footnote>

<footnote>ª See Reimerdes et al 2022 (https://doi.org/10.1088/1741-4326/ac369b) for the TCV Team.</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
Keywords: integrated modelling, ramp-up, QuaLiKiz, TGLF, validation

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

# 1. Introduction

The ramp-up of a Tokamak discharge is the phase where the plasma current is raised from a few tens of Kilo-Amperes to its flattop level. This follows the breakdown and the burn-through, during which the plasma is ionised, and precedes the flattop, where the system reaches the desired plasma current. Understanding and optimizing the ramp is crucial as it influences the duration and characteristics of the flattop. Furthermore, failure in controlling the ramp-up within multiple constraints could result in early termination of the discharge. In particular, plasma control needs to be ensured while remaining within limits such as coil current, coil field, voltage, power and central solenoid force limits.

The control systems benefit from the ability to predict several key plasma quantities before the experiment, among which is worth mentioning the internal inductance $l_i$, the stored energy $W$, the normalized plasma pressure $\beta_N$ and the loop voltage $V_{loop}$. The internal inductance is representative of the current profile peaking and its value needs to be kept low enough (broad current profile) for the coils to ensure plasma position control. For example, the Poloidal Field (PF) system of ITER will only allow a range of $l_i = 0.7–1.0$ [4, 5]. $\beta_N$ and $l_i$, through the shape of the pressure and the plasma current density, are connected to Magneto Hydro Dynamic (MHD) instabilities [6]. Both need to be kept under their stability limits [7, 8]. The stored energy gives an idea on the severity of a possible disruption and $V_{loop}$ is directly related to the flux consumption of the central solenoid.

Impurity content is another crucial factor during the ramp-up phase. It's essential to keep impurity levels, particularly heavy impurities, sufficiently low to prevent excessive plasma cooling from radiative losses. A cooler plasma can trigger stronger resistive modes, which degrade the confinement, and radiate more strongly, cooling the plasma even further. This kind of feedback loops can lead to the loss of plasma.

An ongoing research effort is dedicated to ITER. Experiments aimed at understanding the ITER ramp-up have been performed in various tokamaks. One example is DIII-D, where the impact of various elements on the plasma evolution was examined. These include the impact of the vertical magnetic field $B_z$, the shaping of the plasma during the initial phase, and the potential utilization of Electron Cyclotron (EC) heating early on in the process [9–11].

Predictive modelling of the ramp-up is a great tool in helping to understand the physical mechanisms which determine the evolution of the plasma during this phase. Accurate modelling can in fact not only identify the important processes in each case, but also help in optimizing experiments and control

systems. Two examples are the optimization of the flux consumption and the shape of the current profile at the beginning of the flat-top [12]. Regarding the first, in ITER [13] and future reactors a large fraction of the magnetic flux is consumed during ramp-up and ramp-down, so optimization of these phases can result in longer pulses. On the second, the shape of the current profile has a strong impact on turbulent transport and MHD, and hence on the confinement [14, 15]. The current profile is frozen after the ramp-up by the increase in temperature, for tens of seconds in the planned DEMO. The plasma can then be facilitated in achieving the desired safety factor ($q$) profile, using the current drive systems as actuators, by providing optimized initial conditions for the flattop.

In this paper, integrated modelling of the Tokamak à Configuration Variable (TCV) is presented. TCV is a medium size tokamak, with carbon first wall, maximum magnetic field of $B_T = 1.5$ T, maximum current of $I_p = 1.0$ MA and typical pulse duration around $\tau = 2$ s. The short current diffusion time, especially during the ramp up, allows for $q$ profile relaxation relatively early in the discharge. Distinctive features of TCV are the large flexibility of the magnetic configuration and the very good diagnostic coverage. In particular, the ability to accurately measure the electron temperature $T_e$, the ion temperature $T_i$, the electron density $n_e$ and the carbon density $n_c$ is crucial for the validation exercise that is carried out. Furthermore, the magnetic measurements allow for an estimate of the $q$ profile and of the resistive effective charge $Z_{eff}$. More details on the diagnostics will be presented in section 2.

Modelling of the ramp-up phase is particularly challenging, not only because of the increased nonlinearity that originates from the equilibrium quantities rapidly varying during this phase, but also for the increased difficulty of measurements and for the uncertainty and importance of the initial and boundary conditions. Extensive modelling aimed at understanding the ramp-up was performed during the years, with instances already present during Tokamak Fusion Test Reactor (TFTR) operation [16]. In particular, integrated modelling emerged as one of the most effective tools in studying this transient phase. Different transport models were validated in various machines, enabling validation on a wide range of regimes. Examples can be found in Joint European Torus (JET) modelling with Bohm-GyroBohm (BGB) [4, 17] and the Coppi-Tang (CT) and BGB model for the Experimental Advanced Superconducting Tokamak (EAST) and DIII-D [18, 19].

More recent results include theory-based models such as GLF23 [20], TGLF [2] and QuaLiKiz [1]. The advantage is the higher extrapolability, since a validated theory-based model is not adjusted to the outcome of the experiments in any single or already existing tokamak device, so
it has more general validity and a better founded predictive capability.

Multi-code, multi-machines modelling is particularly effective in providing confidence in extrapolations. Ramp-up modelling was performed at the Axially Symmetric Divertor Experiment (ASDEX) with GLF23 as the quasilinear transport model [21, 22]. JET ramp-up was modelled by leveraging the high speed of the neural network version of QuaLiKiz [23, 24]. Similar modelling was performed at the ‘W Environment in Steady-state Tokamak’ (WEST), focusing on the effect of Nitrogen injection on the plasma evolution [25].

First extrapolations for the ITER ramp-up were performed mostly by scaling the relevant quantities with the Major Radius. Integrated modelling analysis of the ramp-up followed shortly after, focusing a subset of the semi-empirical transport models, such as CT, Neo-Alcator, and BGB [19, 26]. Extrapolation to ITER with integrated modelling and quasilinear models including the ramp-up were also performed [27].

In this work the focus is instead on Low confinements mode (L-modes). This allowed to push the self-consistency and perform, for the first time, integrated modelling with first principle-based-prediction of $j, T_e, T_i, n_e$ and $n_c$ during the ramp up and up to the LCFS. For reasons that will be expanded in section 6, modelling of the TCV ramp-up actually benefits from the self-consistent prediction of the density profiles. To further increase the robustness of this validation, the same settings were used on four different ramp-ups. The platform that was chosen for the modelling is the High-Fidelity Pulse Simulator (HFPS). This python workflow originates from JINTRAC and uses the Integrated Modelling Analysis Suite (IMAS) I/O standard, which supports any tokamak and enables easy comparison with the experimental data.

Such comparison is crucial in identifying important parameters, finding suitable settings and assessing if the agreement between model and experiment is good enough for specific applications. For example, control algorithms for JT-60SA were shown to be able to tolerate perturbations around 10% for $l_{i3}$ and $\beta_N$ [28, 29]. Various metrics can be used for this purpose, each with unique advantages and disadvantages. The wider workflow presented in this paper starts from the TCV experimental data and calculates this agreement, compressing then the information to ensure readability.

Before performing the integrated modelling, however, confirmation is needed that the quasilinear models are able to perform in the specific conditions of the ramp up. From the turbulent transport models point of view, the less validated corner of the parameter space features high $T_e/T_i$, high $q$, high collisionality and high normalized density gradient $R/L_n$. The paper is then structured as follows:

In section 2 the experimental discharges under investigation will be described. In section 3, dedicated gyrokinetic and MHD analysis of the same shot will be presented. The details of the integrated modelling and the respective results are reported in sections 4 and 5. Sensitivities on the numerous choices regarding the modelling inputs and parameters will be presented in section 6. A few metrics of agreement will be

presented and used to assess the impact of these sensitivities in section 7.

## 2. Description of the experimental shots

The four discharges analysed in this paper are all ohmic, inner-wall limited within the modelled time-interval, with deuterium as main ion. Shaping and $Z_{eff}$ are generally low, the flattop current, ramp rates and line average densities vary among the discharges. Shaping and global parameters at the beginning of the modelling and at flattop are listed in table 1.

The HFPS assumes a fully ionized plasma and it is not equipped to evolve a partially ionized plasma. The physics of the burn-through phase [30] is therefore not captured adequately in the modelling. The first High Resolution Thomson Scattering (HRTS) measurement was made at 0.034s for all the experiments, which makes initial conditions available just after the burn-through. For these reasons the simulations were started at 34 ms after breakdown. HRTS measurement were then taken every 17 ms. It was noticed that one of the HRTS laser paths was misaligned. Because of the slightly less correct alignment, the intersection the laser path with the observation lines of sight was deteriorated and therefore the detector received a weaker portion of the scattered light. This lead to an overestimation of the density profile every ~50 ms. This modulation was not present on the Far Infrared interferometer (FIR) measurements, which were then used to renormalize the HRTS profiles. The density at the separatrix from the fit to the HRTS was observed to be occasionally high at the beginning of the discharge in conjunction with bad data quality in that region. For these few points, the value at the separatrix was limited to $10^{19}m^{-3}$.

The Diagnostic Neutral Beam Injector (DNBI) was available for these discharges, allowing measurements of $T_i$ and of the density of the sixth ionization state of carbon $n_C^{6+}$. These measurement are not often available during the ramp-up and contribute to extend the robustness of the validation. The time of the first measurement varies shot by shot, with example of measurements as early as 44 ms after breakdown. Measurements are taken every 24 ms. The three lines of sight of the spectrometers looking at the DNBI in TCV are referred to as SYS01-2-3. For Shot #64965, the active to passive ratio for SYS02 was low and the quality of the data was not high enough, so only data from SYS01 and SYS03 were included in the estimation of $T_i$. Only SYS01 had data with sufficient quality for the carbon density. The measured quantity refers to the last ionization state of Carbon, which is a good approximation of the total carbon density in the inner core, but might be an underestimation towards the separatrix.

All the data were stored using the Interface Data Structure (IDS), within the Integrated Modelling and Analysis Suite (IMAS) workflow [31]. In-house matlab routines were used to extract the TCV data and save them in the IDS structures. The IDSs could then be given directly as input to the HFPS.
<table><caption>**Table 1.** Total $I_p$, elongation $\kappa$, triangularity $\delta$, and effective charge $Z_{eff}$ during the flattop for the studied discharges. The starting time for the flattop phase is reported in the first row, to emphasize the different ramp rates. In the cells for which two values are provided, the first indicates the minimum and the second the maximum value within the considered time interval.</caption><thead><tr><th>Shot #</th><th>64 862</th><th>64 954</th><th>64 958</th><th>64 965</th></tr></thead><tbody><tr><td>t<sub>start</sub> (flattop) (s)</td><td>0.1</td><td>0.15</td><td>0.15</td><td>0.3</td></tr><tr><td>I<sub>p</sub> (MA)</td><td>0.08–0.19</td><td>0.07–0.25</td><td>0.07–0.25</td><td>0.08–0.32</td></tr><tr><td>κ (-)</td><td>1.0–1.5</td><td>1.0–1.5</td><td>1.0–1.5</td><td>1.0–1.5</td></tr><tr><td>δ (-)</td><td>−0.04–0.2</td><td>−0.04–0.2</td><td>0.0–0.2</td><td>−0.04–0.2</td></tr><tr><td>n<sub>e</sub> (line ave) (10<sup>19</sup>m<sup>−3</sup>)</td><td>2.5–2.5</td><td>1.0–3.0</td><td>2.0–4.0</td><td>2.0–5.5</td></tr><tr><td>Volume (m<sup>3</sup>)</td><td>0.6–1.3</td><td>0.6–1.4</td><td>0.5–1.4</td><td>0.7–1.3</td></tr><tr><td>Z<sub>eff</sub> (-)</td><td>1.2</td><td>1.2</td><td>1.2</td><td>1.2</td></tr><tr><td>B<sub>T</sub> (T)</td><td>1.4</td><td>1.4</td><td>1.4</td><td>1.4</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 1.** Traces for the reference shot #64965. Plasma current is reported in (a), $I_{i,3}$ in (b), $v_{loop}$ in (c), $q$ at $\rho_{tor} = 0.95$ in (d) and the electron temperature on axis in (e). Traces for the input and radiated power are shown in (f), in blue for the ohmic heating, in purple for the DNBI, in yellow for the radiation and in green for the total heating power.</figcaption></figure>

The raw data and the errorbars were stored inside the same IDS structures, allowing a certain level of automation of the entire validation framework.

A set of time traces for the reference shot for the modelling, #64965, are shown in figure 1.

### 3. Higher fidelity turbulent transport and MHD modelling

#### 3.1. Turbulent transport modelling

A necessary step to enable integrated modelling is the validation of the quasilinear models against higher fidelity

gyrokinetic in the relevant parameter space. In the ramp-up case this parameter space is characterized by high $T_e/T_i$, high $q$, high collisionality, low plasma shaping, low beta, and relatively low effective charge ($Z_{eff} \sim 1.2$). A linear stability analysis was performed with the flux-tube delta-f gyrokinetic code GKW [32] at two radial positions, $\rho = 0.46, 0.65$ and three different times, $t = 0.068, 0.118, 0.164(s)$. The physics parameters were taken from the experimental fits to discharge #64965 and are reported in table 2. Here $\rho$ is $\rho_{tor}$ and is defined as $\rho = (\phi_{tor}/\phi_{tor,LCFS})^{1/2}$, where $\phi_{tor}$ is the toroidal magnetic flux. $\rho$ will be use in place of $\rho_{tor}$ for the rest of the paper.

The growth-rate spectra are reported in figure 2. Despite the high collisionality, the turbulence is dominated by $R/L_n$
<table><caption>**Table 2.** Physical parameters used for the linear and nonlinear simulations. Values are taken from the fits to the experimental data for shot #64965. Variables without specified dimension are dimensionless. $\nu_{eff}$ is the effective collisionality [33] and $R/L_J$ is the normalized gradient of the quantity $J$.</caption><thead><tr><th></th><th colspan="3">0.46</th><th colspan="3">0.65</th></tr></thead><tbody><tr><td>ρ</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>t (s)</td><td>0.068</td><td>0.116</td><td>0.164</td><td>0.068</td><td>0.116</td><td>0.164</td></tr><tr><td>q</td><td>1.8</td><td>1.3</td><td>1.2</td><td>2.3</td><td>1.8</td><td>1.6</td></tr><tr><td>s</td><td>0.4</td><td>0.6</td><td>0.6</td><td>0.9</td><td>1.2</td><td>1.3</td></tr><tr><td>κ</td><td>1.1</td><td>1.2</td><td>1.3</td><td>1.1</td><td>1.2</td><td>1.3</td></tr><tr><td>β(10<sup>-2</sup>)</td><td>1.5</td><td>3.0</td><td>4.5</td><td>0.6</td><td>1.1</td><td>1.8</td></tr><tr><td>ν<sub>eff</sub></td><td>1.4</td><td>1.3</td><td>1.6</td><td>3.5</td><td>3.4</td><td>3.6</td></tr><tr><td>T<sub>e</sub>/T<sub>i</sub></td><td>1.5</td><td>1.3</td><td>1.2</td><td>2.1</td><td>1.4</td><td>1.4</td></tr><tr><td>R/L<sub>T<sub>i</sub></sub></td><td>6.2</td><td>8.8</td><td>8.3</td><td>5.5</td><td>7.5</td><td>5.5</td></tr><tr><td>R/L<sub>T<sub>e</sub></sub></td><td>11.3</td><td>11.2</td><td>11.7</td><td>9.8</td><td>9.8</td><td>9.2</td></tr><tr><td>R/L<sub>n</sub></td><td>5.8</td><td>7.7</td><td>7.5</td><td>4.9</td><td>5.0</td><td>4.6</td></tr><tr><td>n<sub>C</sub>(%)</td><td>0.8</td><td>0.9</td><td>0.9</td><td>0.7</td><td>0.8</td><td>0.8</td></tr><tr><td>R/L<sub>n<sub>C</sub></sub></td><td>6.9</td><td>7.2</td><td>6.4</td><td>5.2</td><td>0.9</td><td>1.4</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Normalized growth rates $\gamma$ as predicted by GKW as a function of the normalized poloidal wavenumber $k_{\theta}\rho_i$. All calculations are made for shot #64965. Two radial positions are shown, $\rho = 0.46$ ($r/a = 0.5$) on the left plot and $\rho = 0.65$ ($r/a = 0.7$) on the right plot. For each, three time-points at $t = 0.068, 0.118, 0.164(s)$ are plotted in red, blue and purple. Open circles are used when the frequency is in the electron drift direction and crosses when the frequency is in the ion drift direction. The predicted turbulence regime is TEM for $t = 0.068$, transitioning to mixed ITG-TEM for $t = 0.118, 0.164(s)$. ETG is unstable for $\rho = 0.65$ and $t = 0.068(s)$.</figcaption></figure>

driven TEM, especially in the early phase. At higher $k_{\theta}\rho_i$ ETG is predicted to be unstable, but only in the earliest and most outward case. As $s/q$ increases with time, ETG is stabilized and the TEM regime transitions to a mixed TEM-ITG.

Nonlinear gyrokinetic simulations were then performed for the same parameters. Details on the GKW code resolution are reported in table 3. The simulations were electromagnetic, with both $\delta A_{\parallel}$ and $\delta B_{\parallel}$ taken into account. Collisions and flow shear (parallel and perpendicular) were included and the Miller parametrization was chosen for the equilibrium. The fluxes were time averaged over $> 200R/v_{thi}$ after the non-linear overshoot. Inspection of the fluxes spectra reveals that most of the heat flux arises in the range $k_{\theta}\rho_i \sim 0.4-0.6$ which corresponds to the peak in the fluctuation intensity.

The position of this peak is relatively well captured by the usual quasi-linear estimate $\gamma/\langle k_{\perp}^2 \rangle$. To confirm the applicability of quasilinear models in this regime, the non-linear fluxes predicted by GKW were compared to the quasi-linear prediction of TGLF. For this standalone comparison, TGLF simulations were performed with the saturation model SAT2 (NBASIS_MAX = 6, NMODES = 2, NKY = 11) to match the integrated modelling section settings. The TGLF fluxes, normalised to the gyro-Bohm fluxes following IMAS conventions, are shown as a function of the GKW fluxes in figure 3 for the 6 cases described above. Except for a single case for which GKW fluxes are around 4 times larger than TGLF ones, the heat and particle fluxes predicted by TGLF are in the range expected for this kind of exercise, giving confidence the model
<table><caption>**Table 3.** Numerical parameters used for the nonlinear simulations. $k_\theta$ and $k_f$ are the radial and poloidal wave numbers, $\rho_i$ is the ion gyro-radius. NX, NY and NS are the number of grid points in the radial, binormal and parallel direction respectively.</caption><thead><tr><th>NX</th><th>k<sub>f</sub>ρ<sub>i</sub><sup>min</sup></th><th>k<sub>f</sub>ρ<sub>i</sub><sup>max</sup></th><th>NY</th><th>k<sub>θ</sub>ρ<sub>i</sub><sup>min</sup></th><th>k<sub>θ</sub>ρ<sub>i</sub><sup>max</sup></th><th>NS</th></tr></thead><tbody><tr><td>339</td><td>0.04–0.0</td><td>7.3–8.0</td><td>31</td><td>0.05</td><td>1.5</td><td>32</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 3.** Quasi-linear TGLF versus non-linear GKW ion heat (a), electron heat (b), electron particle (c) and carbon particle (d) fluxes for the 6 cases: $\rho = 0.46 (r/a = 0.5)$ and $\rho = 0.65 (r/a = 0.7)$ at $t = 0.068, 0.118, 0.164(s)$.</figcaption></figure>

can be applied in integrated modelling simulations, including the prediction of particle and impurity transport. As often, the fluxes ratio are much better matched, see figure 4 suggesting that the results could be improved by developing a more accurate saturation model.

As shown in figure 4, the electron to ion heat flux ratio is very large at the first time slice, $Q_e/Q_i \sim 6$, when the turbulence is deep in the TEM regime. It then quickly drops to a more moderate value but remains above unity.

## 3.2. MHD modelling

MHD modes are not included in the integrated modelling. When MHD modes are unstable and lead to a sizable modification of the profiles, disagreement between the experiments and the integrated modelling is to be expected. While it is relatively easy to identify MHD unstable modes from the magnetic measurements, predicting them before the experiment is extremely challenging. This is mostly because even small differences in the evolution of the plasma can lead to the destabilization of a mode.

For example, discharges #64954 and #64958 are similar, but #64954 features a locked mode while #64958 does not. The mode significantly redistributes the current, while the effect on confinement is not easy to disentangle from the effect of the initial conditions. The two *q* profiles reconstructed by the equilibrium reconstruction code LIUQE [34] code are very close, but calculations of the classical stability index $\Delta'$ do show that the 2/1 tearing mode is more unstable for #64954. This small difference in the current distribution most likely originates from a small difference in the density and ends up being enough to finally trigger a locked mode at $t \sim 0.6s$. Normalized $\Delta'$ maxima and respective positions are shown in figure 5. The calculations were performed by solving the

tearing equation in cylindrical geometry:

$$ \frac{1}{x} \frac{d}{dx} \left( \frac{xd\psi}{dx} \right) - \left( \frac{m^2}{x^2} + \frac{\mu_0 q (dJ/dx)}{\epsilon B_0 x (1 - nq/m)} \right) \psi = 0 \quad (1) $$

where *x* is the radial coordinate, $\psi$ is the perturbed poloidal flux, $\epsilon$ the inverse aspect ratio, $B_0$ the toroidal magnetic field, *n* and *m* the toroidal and poloidal mode numbers.

Note that, since MHD modes are not included in the integrated modelling, this does worsen the comparisons between model and experiment for these discharges. It would be interesting to use the modelling to predict in advance shots that will or will not develop locked modes, but this requires an accurate knowledge of the initial conditions for the *q* profile. As will be shown in section 6, this is not available yet.

## 4. Input preparation and modelling settings

The same suite of codes described in [35] was exploited to perform the integrated modelling simulations. An important difference is that the HFPS was used instead of JINTRAC. The use of the IMAS I/O standard enables the modelling of any tokamak and was instrumental in this work with TCV. Regarding the settings of the modelling, various decisions had to be taken to avoid numerical instabilities, nonphysical behaviour originating from poor data quality and ensure good agreement with the experimental data.

The simulations were set up to self-consistently predict the current density *j*, the electron temperature $T_e$, the ion temperature $T_i$, the electron density $n_e$ and the carbon density, $n_C$. The neoclassical transport, crucial since it controls the diffusion of the current, was predicted with NCLASS. Collisions and ETG scales were taken into account in QuaLiKiz, which is electrostatic. $B_\perp$ effects were included in TGLF, while $B_\parallel$ effects were turned off. Numerical settings for TGLF are listed
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Comparison of the electron to ion heat flux ratio (a) and particle to ion heat flux ratio (b) predicted by TGLF (blue squares) and GKW (red triangles) as a function of time at $\rho = 0.46$ ($r/a = 0.5$) and $\rho = 0.65$ ($r/a = 0.7$).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 5.** Calculation of stability for the 2/1 tearing mode. The classical stability index $\Delta'$, multiplied by the minor radius $a$, is shown in (a). Discharge #64958, which does not develop a locked mode, is shown in blue, while #64954, which does develop a locked mode, is shown in orange. The position of the $\Delta'$ maximum in $\rho$ is shown in (b) with the same colour-scheme.</figcaption></figure>

<table><caption>**Table 4.** Numerical parameters used for TGLF simulations within integrated modelling.</caption><thead><tr><th># basis funct.</th><th># ETG k<sub>y</sub> modes</th><th># eigenmodes</th><th>sat rules</th><th># grid pts. in JETTO</th></tr></thead><tbody><tr><td>6</td><td>8</td><td>2</td><td>2</td><td>25</td></tr></tbody></table>

in table 4. The toroidal momentum is usually very low in TCV ohmic discharges and was set to zero. For both QuaLiKiz and TGLF the impact of the background $E \times B$ shear was ignored. This effect is important close to the separatrix, but the lack of reduced models for both the precise value of the radial electric field $E_r$ and its effect on the turbulence prevents its self-consistent inclusion.

As stated in section 1, the modelling was initiated 34 ms after breakdown for all the discharges. The transport was predicted using QuaLiKiz [1] for most of the simulations, because of the lower computational cost, but TGLF [36] was deployed for comparison. The QuaLiKiz neural-networks [24, 37] were available, but preliminary runs showed that the TCV parameter space does not lie within the databases used for the training of the neural networks, thus making them not suitable for this application. Specifically,

The magnetic equilibrium reconstruction performed by the LIUQE code and stored in the IDS structure was used to impose the shape of the LCFS. The Equilibrium and Stability COde (ESCO) [38] was then used to self-consistently evolve the equilibrium within the LCFS boundaries. Predicting the experimental current distribution with high accuracy requires the self-consistent modelling of the sawteeth, which avoids too high current density close to the axis. At this end, the sawteeth trigger was calculated with the Porcelli model [39] and assuming a full re-connection process with the Kadomstev model. Note that the sawteeth frequency is very high for both experiment and model, so the central temperature never has enough time to increase and remains relatively flat in the inner core.

Since the Motional Stark Effect (MSE) diagnostic was not available for this discharge, a constrained reconstruction of the q-profile was not available. The error on the initial conditions
for the q profiles was therefore very large. Since the experimental appearance of the first sawtooth was observed shortly after $t=0.034$, the reconstructed q profile was rescaled to impose $q_{axis,t=0.034} = 1.15$. $q_{95,t=0.034}$ was kept constant.

The number of neutrals that cross the separatrix to then be ionized is not measured directly. At the same time, no model capable to calculate it while being computationally cheap enough to be coupled in a set of dynamic integrated modelling simulation is available yet. For these reasons, no value for the number of particle crossing the separatrix is available. To face this issue, a feedback control was imposed to set the neutral source through a close feedback loop on the line averaged density measured by the interferometer. This is analogous to what is done in the experiments and was already demonstrated previously for ASDEX L mode modelling [40]. The ionization energy per atom was set at 13.6 eV, while the neutral energy was set to 30 eV. The latter is assumed to be the energy with which the neutrals cross the separatrix and influences the neutral penetration. Since the recycling was set to zero, the feedback is performed directly over the total number of particles crossing the separatrix.

Due to its prevalence, carbon was chosen as the only impurity. Other impurities might be present, but not in concentrations large enough to significantly modify dilution and $Z_{eff}$. A boundary for the carbon density was set at $\rho=1$. An internal boundary condition fixing the carbon density inside the separatrix can result in numerical problems that can significantly influence the evolution of the early ramp-up when the initial impurity content is large. This is not the case for the reference simulations, but can be the case for the sensitivities.

Instead, the boundary conditions were set by setting an escape velocity of $2 \text{ m s}^{-1}$, zero recycling and an effective carbon puff. The carbon incoming flux was normalized on a shot-by-shot basis using the $Z_{eff}$ measured during the flattop phase and the line-averaged density specified for the feedback loop. Specifically, the carbon puff $\Gamma_C$ was calculated as

$$ \Gamma_C = \Gamma_C(\text{ref}) \frac{n_{e,ave}}{n_{e,ave}(\text{ref})} (Z_{eff} - 1) \quad (2) $$

where $\Gamma_C(\text{ref}) = 4 \times 10^{18} \text{m}^{-3} \text{s}^{-1}$ and $n_{e,ave}(\text{ref}) = 10^{19} \text{m}^{-3}$. In the future it could be interesting to improve this simple scaling or to modify the impurity influx in time to match the whole trace of the $Z_{eff}$ or carbon measurements. The initial $n_C$ profiles were set to impose a $Z_{eff} = 1.5$. This value is slightly higher than the flattop values to avoid assuming a large impurity influx at the start of the simulation, which would be needed to maintain a constant $Z_{eff}$ with an increasing $n_e$.

For the $T_e$, $T_i$ and $n_e$ the boundary was defined at the LCFS. A number of choices had to be made for the initial and boundary conditions when the experimental data were not available or not of sufficient quality. As mentioned above, the initial conditions for the safety factor profile were chosen to impose $q=1.15$ at the magnetic axis. $Z_{eff}$ at $t=0.034(\text{s})$ value was chosen as 1.5. To ensure the physical relevance of the results and filter fits to low quality or sparse experimental data, limits were set for $T_e$ and $n_e$ at the separatrix. These limits were set to $25(\text{eV}) < T_{e,sep} < 150(\text{eV})$ and $n_{e,sep} < 10^{19} \text{m}^{-3}$. The

density ceiling was only surpassed at the beginning of one discharge in conjunction with poor experimental data coverage and should not be interpreted as a hard limit. Furthermore, the $T_i/T_e$ ratio was restricted to 3.0, which is closer to SOLPS modelling results [41], even in cases in which measurements at the boundaries exceeded this value.

Since the simulations start at the same time as the first HRTS measurement, the initial $T_i$ profile was chosen to be the same as the first actually measured $T_i$ profile. This might be a strong assumption in cases when the first available measurement is taken ~100(ms) later, but the sensitivity to the initial conditions for the $T_i$ profile was observed to be modest on the evolution of $T_i$ itself and negligible on $T_e$ and $n_e$. Due to low quality of the $T_e$ measurements close to the separatrix for the first few HRTS data acquisition, $T_e$ is imposed at the separatrix at $t=0$ and then interpolated linearly to its measured value at $t=0.1(\text{s})$. $T_{e,t=0}$ is chosen arbitrarily to be 60(eV), but sensitivities on this value are performed in section 5. Note that lower values cause the simulations to suffer numerical instabilities and crash for some of the discharges.

Numerical considerations were taken into account to facilitate accurate simulations. A smaller maximum timestep of $10^{-6} \text{s}$ was implemented at the beginning to prevent profile collapse, which can occasionally occur when impurities are present and the FRANTIC code is used. This maximum timestep is then quickly raised to the nominal value of $10^{-3} \text{s}$. Note that the energy confinement for these shots is ~20 ms and that an adaptive timestep scheme is used in JETTO. To avoid fluxes being identically zero, a negligible transport component scaled as Bohm diffusivity was included. Because of numerical stability, the impurity puff was initiated 40 ms after the beginning of the simulation. Additionally, the input IDS slice was updated with a time interval of $10^{-3} \text{s}$ to ensure a smooth behaviour of the equilibrium quantities.

## 5. High fidelity pulse simulator results

Results for the modelling of shot #64965 are shown in figures 6 and 7. Simultaneously comparing multiple channels provides confidence that the self-consistent modelling is effectively capturing the feedback loops between the various quantities. Radial profiles are shown for $T_e$, $T_i$, $n_e$ and $n_C^{6+}$ and refer to $t=0.167 \text{ s}$, which is the latest instant used for the gyrokinetic analysis in section 3. Time traces are shown for the internal inductance $l_i$ and the reconstructed stored energy $W_{tot}$. The modelled profiles are shown on top of the experimental data, raw for the radial profile and reconstructed from the equilibrium reconstruction code LIUQE for the time traces.

Good agreement for this type of modelling is obtained for $T_e$ and $n_e$, with the prediction and the data showing only a ~15% and ~10% difference respectively. The precise method used to obtain these percentages will be explained below. It is somewhat surprising that both QuaLiKiz and TGLF show such agreement with the experiment, particularly further out towards the separatrix.
<figure><img src="image_15.png" /><figcaption>**Figure 6.** Profiles for discharge #64965. In (a) $T_e$ as measured from the High Resolution Thomson Scattering (HRTS) and as predicted by the integrated modelling with the HFPS are shown in blue and yellow respectively. The same is shown in (b) for $n_e$. In (c) and (d) the experimental data for $T_i$ and $n_C^{6+}$ are instead obtained through the CXRS. All profiles are shown for $t = 0.167$ s.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 7.** Time traces for discharge #64965. The internal inductance $l_{i,3}$, and the stored energy $W_{thermal}$ are shown in (a) and (b). The values as reconstructed by LIUQE and as predicted by the integrated modelling with the QuaLiKiz and TGLF as the turbulent transport model are shown in blue, orange and green respectively.</figcaption></figure>
The resistive modes which are in fact expected in this region [42] are not included in QuaLiKiz and excited only at very large collisionality in TGLF. We in fact do observe instances of overestimation of the temperature gradients close to the LCFS, especially in QuaLiKiz. In most of these cases, however, the mismatch in gradient is limited to a region narrow enough that the difference in temperatures and density remains limited. Discharges with higher density or, especially for QuaLiKiz, more extreme shaping might show higher disagreement.

The profile for $T_i$ seems instead to be underestimated. Closer to the boundary conditions this effect might be due to $E \times B$ physics, which is ignored here. If this results in an overestimation of the particle transport, the approximation could also increase the charge-exchange ion heat sink calculated by FRANTIC, further decreasing the ion temperature. A second explanation is that ignoring the DNBI is too strong of an assumption and the modest heating, around 75 kW in this case, is enough to raise $T_i$ by the missing ~50 eV. Modelling with NBI is however quite complex and will be left for future work. For $n_C^{6+}$ the gradients are not very well captured, especially close to the LCFS. Note that the total carbon density is not in feedback and the settings are the same for both QuaLiKiz and TGLF, so some difference in total concentration is expected.

The mismatch of $l_i$ for the earliest times might be due to the choice of the initial conditions. The sensitivity of this and other quantities to the initial $q$ profile will be studied in section 6. The energy content is slightly underestimated for QuaLiKiz and almost perfect with TGLF, even thought it is obtained with an overestimated $T_i$ and an underestimated $T_e$. The experimental value of $W_{thermal}$ at $t=0.034$ s is too low for a precise value to be provided by the LIUQE reconstruction. For this first point the energy can be calculated from the measured $T_e$, $T_i$ and $n_e$ and it does compare quite well with the modelling. However, since the values of $W_{thermal}$ within the integrated modelling are calculated as the integral of $T$ and $n$, using $W_{thermal}$ as an accuracy estimator does not add much information with respect to the comparison of the profiles. This is the reason behind the choice to show the energy content calculated by LIUQE in figure 7(b).

In figure 6 the agreement is shown for one specific timeslice. A strategy was then devised to illustrate the quality of the prediction through the whole ramp-up. A simple measure of agreement defined as

$$ \sigma_t = \sum_{i=0}^{N} \frac{2}{N} \left| \frac{d_{fit, \rho_i} - d_{model, \rho_i}}{d_{fit, \rho_i} + d_{model, \rho_i}} \right| \quad (3) $$

was calculated. Here, $\rho_i$ refers to one radial position, while $N$ represents the total number of radial positions considered. $d_{fit}$ indicates the fit to the experimental data, while $d_{model}$ indicates the modelling data extrapolated at the time and radial position of the corresponding point of the fit.

The measure is of simple interpretation and describes the average relative distance between the modelling and the fit to the experimental data. When multiplied by 100, [3] can be interpreted as an average percentage agreement. This choice is arbitrary and other measures of (dis)agreement will be presented in section 7.

Time traces of this measure are shown in figure 8, Facilitating a rapid assessment of regions of good and bad performance of the models. The prediction for $T_e$ and $n_e$ is in general of good quality, particularly after ~100 ms. The disagreement before this time might be due to three reasons. The first might come from the fact that the disagreement lasts a few confinement, and current diffusion, times. This points in the direction of the initial conditions being important. In particular, the initial current profile is imposed and not measured. The real initial current distribution is therefore not known precisely. A second reason could be that QuaLiKiz and TGLF are not able to capture the transport this early in the discharge, with low $T_e$ and low volume. Lastly, MHD modes not calculated by the model could be unstable.

The prediction of $T_i$ is instead best at the start of the simulation, but progressively deteriorates up to a ~30% disagreement with the fit to the experimental data. Note that $T_i$ was limited to $3T_e$ and part of the disagreement originates from the boundaries. The boundaries have a large weight in formula 3 due to the small denominator, so this contribution can be relatively large. For shot #64965 the removal of this limit brings the disagreement down to ~5%. Interestingly, this has a very limited impact on the $T_i$ value on axis. These simulations, although improved in agreement, rely on physically unrealistic values of $T_i/T_e$ and instead suggest an inconsistency in the measurements.

The slight underestimation of the total stored energy $W_{tot}$ follows from the underestimation of $T_i$. $l_i$ is generally overestimated, especially by TGLF, but the prediction is still accurate for 10%. As mentioned in section 4, $V_{loop}$ is underestimated for the first ~100(ms), probably due to the presence of non-carbon impurities. No attempt was made to match the experimental radiation. The radiated power depends on the total carbon content, which in this modelling is only roughly matched through the constant carbon inflow. Furthermore, proper radiation modelling would need the inclusion of heavier impurities, which increases the complexity of the simulation and its left for future work. Similarly, the heating power and particle source originating from the DNBI is ignored.

The power balance shows $Q_e/Q_i \sim 3$ at the start of the simulation, slowly falling to ~1 as the density rises and the power exchange between electrons and ions increases. It is interesting to note that a non-negligible fraction of the ion heat transport is due to the neoclassical transport. Outside the sawtoothing region, $\chi_{nc} \sim 0.5(m^2s^{-1})$. Regarding the turbulence regime, QuaLiKiz, TGLF and GKW qualitatively agree. The plasma is TEM unstable across the whole radius. For $\rho \sim < 0.5$ ITG is stable at the beginning of the simulation and then gradually destabilizes. Moving to inner radii, ITG is gradually stabilized and the plasma remains in a pure TEM regime for the whole duration of the modelling. Qualitatively, this agrees with the higher fidelity modelling performed in section 3.

Lastly, a more direct comparison with the experimental current profile was achieved by comparing the time of the first sawteeth as well as the radial position of the inversion radius. On the experimental side, the inversion radius was calculated from the soft x-rays emission, while for the modelling the inversion radius was taken as the point where $q = 1.1$.
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Time trace of the measure of agreement described in formula 3. The evolution of $\sigma_t$ for $T_e$, $n_e$, $T_i$ and $l_{i3}$ is shown in (a), (b), (c) and (d) respectively. The yellow lines show the measure calculated on HFPS-QuaLiKiz, while a comparison with HFPS-TGLF is represented by the green lines. The black lines indicate the time for which the profiles in figure 6 are shown.</figcaption></figure>

The modelling error bars were chosen to cover the interval $1.05 < q < 1.15$, so that the large confidence interval ensures a fair comparison. The results for the reference simulation of shot #64965 and two sensitivities are shown in figure 9. For sensitivities, the initial $q$ profile was modified to keep $q_{95}$ constant and impose $q_{axis,t=0.034} = 1.6, 1.0$. The evolution is similar in all cases and agrees well with the experimental data. A delay of ~ -10ms, ~20 ms and ~40 ms is observed, for $q_{axis,t=0.034} = 1.0, 1.15$ and 1.6 respectively.

**6. Integrated modelling sensitivities**

After assessing the validity of the modelling on one TCV ramp-up, the relevance of this study was enhanced by running the sensitivity on all four discharges described in table 1. An extensive set of simulations was then run to assess the statistical sensitivity of the modelling to physical and numerical settings. Each time, one specific aspect was modified with respect to the reference simulation. Four sensitivities are reported in

figure 10 for shot #64862, for which some of the trends are clearer. A different set is reported in figure 11 using the metric described in section 5. For this scope, the metric was averaged over time and condensed to one single number for each quantity. For both figures, the relative simulations will be highlighted explicitly when describing them. Not all sensitivities mentioned in the text are shown in figures, particularly when the impact is low and the profiles are very close to the reference profiles.

Such plots can be used to explore the impact of the choices for the modelling settings, study the physics and eventually to further refine the prediction. These sensitivities are TCV specific and other machines might behave slightly differently, but the validation for high $q$, high collisionality and high $T_e/T_i$ is transferable to other tokamaks. The main differences when modelling different machines are the timescale for the current diffusion, longer in larger tokamaks, and the impurity composition, which can have a significant impact when the tungsten concentration is large. This exercise should then be repeated for a larger, tungsten machine. This would however
<figure><img src="image_18.png" /><figcaption>**Figure 9.** Comparison of experimental (blue) and modelled (green) inversion radius for shot #64965. The reference simulation for discharge #64965 is shown in (a). Two sensitivities where the *q* profile was manipulated to have a different initial *q*<sub>axis</sub> are shown in (b) and (c). The calculation of the values and error bars is described at the end of section 5. The best agreement is obtained for the run with the lowest initial *q*<sub>ρ=0</sub>.</figcaption></figure>

be simplified by the use of IMAS, which ensures the re-
usability of the software developed in this work.

**6.1. Boundary conditions**

The boundary conditions were modified by increasing or
decreasing the value at the boundaries for T<sub>e</sub>, T<sub>i</sub> and n<sub>e</sub> and
then interpolating the profile to have the same value on the
axis. This does not have a huge effect on the plasma state. The
highest sensitivity is of T<sub>e,sep</sub> on l<sub>i3</sub>, up to 5% when increasing
T<sub>e,sep</sub> by 40 eV, which corresponds to its value at the top of
the error bar. This change is indeed in the expected direction,
with a higher temperature at the separatrix leading to a lower
resistivity in that region and therefore to a flatter q profile. T<sub>i</sub>
has an even more limited impact on the profiles. This impact is
mainly due to the coupling with FRANTIC. The temperature
of the neutrals is kept constant, so more energy is required to
heat the neutrals to this higher temperature and the ion heat
sink increases. As a result, the ion temperature at ρ = 0.9 stays

roughly constant independently from its value at the separat-
rix. The feedback loop on the density leads to a similar reas-
oning on n_e, at least for variations below ~30%. Interestingly,
R/L_n_e and R/L_n_c close to the boundary conditions seem to con-
sistently increase with lower T_i/T_e.

6.2. *Initial conditions*

Initial conditions for temperature and density have a limited
impact. The profiles converge to the reference simulation 10-
20 (ms) after the difference in the input disappears. The initial
$q$ profile has instead a fairly strong impact on the evolution.
Figure 12 shows the impact of $q_{t=0.034}$ on $l_{i3}$ and highlights
how most of the early disagreement in fact originates from the
choice of the initial conditions. However, the optimal initial
conditions vary from discharge to discharge. Furthermore, the
simulation with the best agreement for $l_i$ does not have the best
agreement with the experimental sawteeth timing. A probable
<figure><img src="image_19.png" /><figcaption>**Figure 10.** A subset of the sensitivities described in section 6, all based on the integrated modelling with the HFPS of shot #64862. $T_e$, $n_e$, $T_i$ and $n_C^{+6}$ are shown in (a) (b) (c) and (d) respectively. The reference simulation for #64862 is shown in blue. A simulation with the electron density imposed from the experiment instead of being predicted is reported as 'Interpretive $n_e$' in yellow. For the simulation in green $P_{rad}$ was set to 200 kW and the impact of ETG was removed in the simulation shown in red.</figcaption></figure>

reason is that the shape of the *q* profile is important and not optimized in these simulations.

Lastly, $Z_{eff} = 4$ at $t = 0.034$ was set as an upper limit value at the start of the simulations to try to assess the impact of the possible presence of non-carbon impurities on the modelling. Note that the carbon content is still set mostly by the carbon puff, so $Z_{eff} = 4$ relaxes to a value of 1.5 after ~100 ms. $l_i$ is marginally lower and $v_{loop}$ is only ~0.5 higher. This is mostly due to an increase of ~150 eV in $T_e$, which decreases the resistivity. The impact on $T_i$ and $n_e$ is marginal.

## 6.3. Physics settings

The impact of various physics settings and mechanisms was analysed next, starting with the neutral particle source. Unsurprisingly, when turning FRANTIC off and neglecting the neutral source the density gradient decreases drastically, which leads to a large increase of the electron temperature.

Note that the density never goes to zero since a non zero value is imposed at the separatrix. Increasing the energy of the neutrals crossing the separatrix in FRANTIC to 150eV increases significantly the penetration of the neutrals and switches the sign of the heat sink due to the Charge exchange and the ionization of the neutrals. The reason is that the neutrals become hotter than the plasma, which should not be seen as physical but as an extremal situation. This leads to a lower $R/L_{n_e}$ and a ~20 eV higher $T_e$ close to the separatrix. The increase in $T_i$ is similarly modest, aside from when the heat sink causes a negative $R/L_{T_i}$, which is the case at the beginning of some of the reference simulations. The difference in $T_e$ reduces $l_i$ by a few percent to 20%, depending on the discharge.

Radiated power, even when increased to 60%–80% of the total ohmic power, proved to be mostly inconsequential for $T_e$ and, interestingly, to increase $T_i$. In this case the radiation was set to be constant through the profile and the low impact on $T_e$ is due to stiffness. The impact on $T_i$ seems instead to be due to
<figure><img src="image_20.png" /><figcaption>**Figure 11.** Sensitivities of the time average of the measure of agreement described in formula (3) for $T_e$, $n_e$, $T_i$ and $l_{i3}$ for (a), (b), (c) and (d) respectively. Different colours represent different discharges, while the different sensitivities are described in the x axis. For (b), the measure is non zero for ‘Interpretive $n_e$’ because of interpolation differences in the fit and the model output times and the limit in separatrix density explained in section 2.</figcaption></figure>

non-diagonal transport and to the difference in $T_i/T_e$ close to the boundary conditions.

The Electron Temperature Gradient (ETG) driven turbulence was turned off in QuaLiKiz to assess its impact. The effect depends on the discharge, but in general $T_e$ and $T_i$ are increased. For the most sensitive discharges, including \#64965, the increase can be ~15% on axis. While ETG is predicted to be unstable close to the boundary conditions, the flux in that region is not predicted to be large and most of the difference appears around mid radius. The profiles for this sensitivity on shot \#64862 are shown in figure 10. A greater impact on $T_e$ is found to originate from $Z_{eff}$, controlled through the carbon puff. On axis, the difference can be up to 200(eV) for a $\Delta Z_{eff}$ of 0.7. Turning on the effect of $E \times B$ shear in QuaLiKiz has the desired effect of increasing the ion temperature, but at the same time leads to an overestimation of $R/L_{Te}$ and $R/L_{Te}$, mostly localized close to the boundaries.

The investigation of possible mismatches in the calculation of the transport close to the boundaries was performed by imposing extra heat electron diffusion. The additional transport was imposed as a Gaussian centred at the separatrix. The amplitude was set as $8(m^2s^{-1})$ and the width as 0.1 in $\rho$ coordinates. This approach is similar to what was done in [8]. The agreement measure is referenced as ‘Extra transport’ in figure 11. Although $T_e$ is indeed closer to the experiment, the density peaking becomes larger in the same region and the overall agreement decreases.

The sensitivity involving Sawteeth justifies the extra computational cost derived from the sawteeth model. In the simulation referenced as ‘No Sawteeth’ in figure 6, the Sawteeth model was turned off. $T_e$ and $T_i$ rise in the core, but the difference in current distribution decreases the shear and both temperatures are lower outside mid-radius. The temperature difference at $\rho = 0.5$ is ~15%–20%.
<figure><img src="image_21.png" /><figcaption>**Figure 12.** Time evolution of $l_i$ with different initial conditions for the $q$ profile. $q_{0, axis}$ refers to the value of $q$ at $t = 0.034$ and at $\rho = 0$. $q_{0, axis} = 4.5$ corresponds to a flat starting $q$ profile, while $q_{0, axis} = 1.0$ immediately starts sawtoothing. Those can be considered as upper and lower limits for $q_{0, axis}$ respectively. The experimental measurement is the output $l_i$ of the LIUQE reconstruction.</figcaption></figure>

A second crucial justification for this modelling is provided by the simulation referenced as ‘Interpretive $n_e$’ in figures 11 and 10. The electron density was imposed instead of predicted. The effect of this setting is more subtle. While for some of the discharges the difference in the predicted $T_e$ is only ~5%, for #64958 and #64862 it is above 20%. This is observed concomitant to regions with high $R/L_{n_e}$ and density gradient-driven TEM turbulence. In these cases, a relatively small decrease of $R/L_{n_e}$ would be enough to decisively reduce the transport. When artificially keeping $n_e$ fixed, however, the system is forced to greatly decrease $R/L_{T_e}$, flattening the $T_e$ profile. The agreement with the experiment is lost.

### 6.4. Considerations towards predict first

A last sensitivity was performed in order to increase the predictive capabilities of this modelling and is referenced as ‘Imposed BC’ in figure 6. The boundary conditions for $n_e$, $T_e$ and $T_i$ were imposed as $7 \times 10^{18}$ (particlesm$^{-3}$), 40 eV and 100 eV instead of being taken from the experiment. No attempt was made in optimizing the precise values, which are arbitrary and simply characteristic of typical TCV experimental conditions. The results are not very different from the reference simulations. Most of the disagreement in $n_e$ and $T_e$ is due to the difference close to the boundaries, which is heavily weighted in the measure of agreement 3.

It is noteworthy that the only experimental data used in these last simulations are the line averaged density $n_e(l)$, $I_p$, the equilibrium boundaries and the flattop $Z_{eff}$. In the experiment $I_p$ is imposed, $n_e(l)$ can be controlled in feedback and $Z_{eff}$ is fairly low and similar between different ohmic shots. Given an already modelled equilibrium, therefore, the pipeline presented in this work should be able to reasonably perform predict first simulations of ohmic TCV discharges.

## 7. Alternative measures of agreement

Alternative measures of (dis)agreement were proposed and are presented here. Each new metric provides useful new information and can be tailored for a specific application. Three examples are applied to the same set of sensitivities reported in figure 11. Such procedure aids in the identification of the most relevant assumptions in the modelling, which in turns directs the direction of future efforts.

### 7.1. Comparison with the fit

Firstly, it needs to be noted that the grid that is used to perform the sum in equation (3) is not necessarily equidistant in $\rho$. The default grid is the grid that was used to perform the fit to the experimental data. The choice of allowing for different grids was made to increase the flexibility of the metrics. Most of the times, for TCV, the grid is more dense towards the boundaries. The effect in the formula is that the boundaries are weighted more.

The first alternative measure is then to simply use a regular grid in equation (3). This is shown for the electron temperature in figure 13(a). With respect to figure 11(a), a few minor differences are visible. First, the marginally higher disagreement for TGLF with respect to QuaLiKiz vanishes. The simulation for which plausible values were imposed as boundary conditions shows agreement comparable to the reference runs. This, unsurprisingly, suggest that most of the difference originates from the region close to the boundaries. Lastly, the disagreement is visibly larger when imposing the density profile, sign that the mismatch originates from the core.

A second measure can be defined as

$$ \sigma_t = \sum_{i=0}^{N} \frac{1}{V_i - V_{i-1}} \frac{2 d_{\text{fit}, \rho_i} - d_{\text{model}, \rho_i}}{N d_{\text{fit}, \rho_i} + d_{\text{model}, \rho_i}} \quad (4) $$
<figure><img src="image_22.png" /><figcaption>Figure 13. Sensitivities for $T_e$ with different methods to calculate the disagreement. Different colours represent different discharges, while the different sensitivities are described in the x axis. A regular grid in $\rho$ is used in (a), the errors are weighted with the volume in (b), the difference is calculated without the absolute value in the summation in (c) and the model is compared directly to the experimental data in (d).</figcaption></figure>

where V is the volume associated with the radial point at which the fit was calculated. The results are shown in figure 13(b) and are similar to the reference. This measure is probably more relevant than the one chosen as reference, but it is of less immediate interpretation.

## 7.2. Comparison with the raw experimental data

A metric with a different purpose is obtained by keeping the sign of the differences instead of calculating the absolute value. This gives a quick indication of systematic overestimations or underestimations of the data. It is important to highlight that the two can balance each other and small disagreement here might not mean a good match with the fit to the experimental data. From figure 13(c) the utility of this measure is however evident. QuaLiKiz seems generally overestimate the profiles, while this is not the case for TGLF. The extra transport seems to be too strong for shots #64965 and #64954, for which the temperature decreases to be lower than

the fits. A more precise estimation of the missing transport close to the separatrix can then be made. The overestimation for the imposed boundary condition might simply mean that the value of $T_{e,\rho=1}$ needs to be lowered. Imposing the density causes the electron temperature to drop, consistently with the explanation of $R/L_{T_e}$ having to flatten to balance the drive when $R/L_{n_e}$ becomes too large.

Lastly, the model can be compared directly with the experimental raw data instead that with the fit. The formula then becomes

$$ \sigma_t = \sqrt{\sum_{i=0}^{N} \frac{1}{N} \frac{(d_{\text{exp},\rho_i} - d_{\text{model},\rho_i})^2}{\sigma_{\text{exp},\rho_i}^2}} \quad (5) $$

Here, N represents the number of radial positions at which the measurements were taken, while $d_{\text{exp}}$ indicates the experimental data. $d_{\text{model}}$ indicates the modelling data extrapolated
<figure><img src="image_23.png" /><figcaption>**Figure 14.** Relative distance as calculated from equation (5). The metric calculated by comparing experimental data with the reference simulation is reported in blue, while the comparison with the fit to the experimental data is reported in orange. When the measure is below 1.0 it means that the simulation (or the fit) is within the experimental error-bars. Scattered experimental data can cause the distance to be greater than one even for the fit.</figcaption></figure>

at the time and radial position of the corresponding experi-
mental point. For this measure, the simulation without saw-
teeth and the one with imposed density profile are penalized,
but it is encouraging to note that for both QuaLiKiz and TGLF
the model is basically within the experimental error-bars. This
measure of easy interpretation, but it does have a problem
when the experimental data are of low quality or scattered.
An example is shown in figure 14 for the ion temperature.
The measure is calculated for both the model and the fit to
the experimental data. It is immediately evident that the dis-
agreement is large even for the fit, thanks to the scatter of the
data and the small error-bars of some of the data point. As can
be seen from the correlation between the agreement with the
fit and with the model, equation (5) contains not only the per-
formance of the model, but also the quality of the experimental
data.

More alternatives could be envisioned, mixing the ones
presented above or introducing new concepts. None of the
measures is perfect and the each time the choice might be a
combination of them and will depend on the specific problem
being solved.

8. Conclusions

Four TCV ramp-up phases of ohmic discharges have been
successfully reproduced with integrated modelling up to the
LCFS, with self consistent prediction of *j*, *T*<sub>e</sub>, *T*<sub>i</sub>, *n*<sub>e</sub> and *n*<sub>C</sub>
Despite the complexity, the modelling showed good agree-
ment with the experimental data. The electron and ion densit-
ies were predicted, increasing self-consistency and accuracy at
the same time. The simultaneous comparison of multiple chan-
nels allowed for a more stringent validation, helping against
overfitting of single quantities. The use of identical settings

across all of the simulation shows the robustness of the setup
and provides a powerful workflow for validation.

Broad agreement was found between the turbulence pre-
dicted by the quasilinear model TGLF and the higher fidel-
ity model GKW. The plasma was found to be dominated by
R/L<sub>n<sub>e</sub></sub> driven TEM, with a later transition to a mixed ITG-TEM
regime. Power balance analysis on both integrated and higher
fidelity modelling found Q<sub>e</sub>/Q<sub>i</sub> > 1. Some measures of agree-
ment were proposed to identify regions of poor performance of
the modelling and condense the information in a small number
of plots.

The extensive sensitivity tests guided the identification
of actuators for the simulations. While boundary conditions
proved to have a low impact, excluding large variations in
*n*<sub>e,sep</sub>, initial conditions for the *q* profile and *Z*<sub>eff</sub> were found
to be important for the first ~100 ms of the ramp-up. The cru-
cial role of the density prediction, including the neutral source,
the sawteeth model and of the transport close to the separatrix
was investigated and confirmed the settings choice.

Future work might develop in various directions. First,
some remaining disagreement on the first 100 ms of the com-
parison requires a wider and more systematic study on the
initial conditions, particularly for the current profile. Proper
boundary conditions for *E*<sub>r</sub>, combined with a deeper under-
standing of its impact on the turbulence close to the LCFS,
should allow its inclusion in the modelling.

On the longer term, the same pipeline presented here can be
applied to an arbitrary number of discharges, eventually from
multiple machines, to perform a more comprehensive valida-
tion. This effort would be useful in identifying which assump-
tions and simplifications are too strong and which experi-
mental conditions are challenging for the model. Future work
can then be guided and directed to improve the models when
accuracy is insufficient.
Finally, the framework is ready for a first attempt at predict first of ohmic, L-mode TCV discharges. On the experimental side, the discharge should not terminate due to MHD, which is not included in the modelling, and the control systems should ensure that the references for $I_p$, $n_e$,lineaveraged and shape are closely followed. On the modelling side, the guesses for $n_e$,sep and the initial conditions should be within 30% of the actual experimental value. Good agreement can be expected, provided that these conditions are fulfilled.

## Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200-EUROfusion). The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). The gyrokinetic simulations were performed on the EUROfusion High Performance Computer Marconi-Fusion under the TSVV11 Project WP2D1. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. This work was supported in part by the Swiss National Science Foundation and was carried out in collaboration with the TSVV-11 team.

## ORCID iDs

M. Marin https://orcid.org/0000-0002-7483-4118

Y. Camenen https://orcid.org/0000-0003-1361-3722

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

F.J. Casson https://orcid.org/0000-0001-5371-5876

R. Coosemans https://orcid.org/0000-0001-8110-3156

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

P. Maget https://orcid.org/0000-0002-8746-8102

A. Najlaoui https://orcid.org/0009-0006-7169-841X

O. Sauter https://orcid.org/0000-0002-0099-6675

## References

[1] 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

[2] 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

[3] Tamor S. 1981 ANTIC: a code for calculation of neutral transport in cylindrical plasmas *J. Comput. Phys.* **40** 104119

[4] Hogeweij G.M. et al 2014 ITER-like current ramps in JET with ILW: experiments, modelling and consequences for ITER *Nucl. Fusion* **55** 013009

[5] Sips A.C. et al 2009 Experimental studies of ITER demonstration discharges *Nucl. Fusion* **49** 085015

[6] Buratti P. et al 1997 MHD activity in FTU plasmas with reversed magnetic shear *Plasma Phys. Control. Fusion* **39** B383

[7] Gribov Y., Humphreys D., Kajiwara K., Lazarus E.A., Lister J.B., Ozeki T., Portone A., Shimada M., Sips A.C. and Wesley J.C. 2007 Progress in the ITER physics basis, chapter 8: plasma operation and control *Nucl. Fusion* **47** S385

[8] Imbeaux F. et al 2011 Current ramps in tokamaks: from present experiments to ITER scenarios *Nucl. Fusion* **51** 083026

[9] Jackson G.L. et al 2008 ITER startup studies in the DIII-D tokamak *Nucl. Fusion* **48** 125002

[10] Jackson G.L., Casper T.A., Luce T.C., Humphreys D.A., Ferron J.R., Hyatt A.W., Leuer J.A., Petrie T.W., Turco F. and West W.P. 2009 Simulating ITER plasma startup and rampdown scenarios in the DIII-D tokamak *Nucl. Fusion* **49** 115027

[11] Jackson G.L., Politzer P.A., Humphreys D.A., Casper T.A., Hyatt A.W., Leuer J.A., Lohr J., Luce T.C., Van Zeeland M.A. and Yu J.H. 2010 Understanding and predicting the dynamics of tokamak discharges during startup and rampdown *Phys. Plasmas* **17** 056116

[12] Felici F. and Sauter O. 2012 Non-linear model-based optimization of actuator trajectories for tokamak plasma profile control *Plasma Phys. Control. Fusion* **54** 025002

[13] Luce T.C., Humphreys D.A., Jackson G.L. and Solomon W.M. 2014 Inductive flux usage and its optimization in tokamak operation *Nucl. Fusion* **54** 093005

[14] Hoang G.T. et al 1998 Magnetic shear effects on confinement and electron heat transport in torre supra discharges with electron heating *Nucl. Fusion* **38** 117–32

[15] Citrin J. et al 2012 Predictive analysis of q-profile influence on transport in JET and ASDEX Upgrade hybrid scenarios *Plasma Phys. Control. Fusion* **54** 065008

[16] Jardin S.C., Bell M.G. and Pomphrey N. 1993 TSC simulation of Ohmic discharges in TFTR *Nucl. Fusion* **33** 371–82

[17] Nunes I., Saibcne G., Lomas P.J., Sips A.C., Lonroth J. and Voitsekovitch I. 2011 ITER similarity current ramp-down experiments at JET *38th EPS Conf. on Plasma Physics 2011, EPS 2011 - Europhysics Conf. Abstracts* vol 35 pp 1328–31

[18] Liu L., Guo Y., Chan V., Mao S., Wang Y., Pan C., Luo Z., Zhao H. and Ye M. 2017 Validating the energy transport modeling of the DIII-D and EAST ramp up experiments using TSC *Nucl. Fusion* **57** 066033

[19] Casper T.A., Meyer W.H., Jackson G.L., Luce T.C., Hyatt A.W., Humphreys D.A. and Turco F. 2011 Validation of the thermal transport model used for ITER startup scenario predictions with DIII-D experimental data *Nucl. Fusion* **51** 013001

[20] Waltz R.E. et al 2013 A gyro-Landau-fluid transport model *Phys. Plasmas* **24** 2482–96

[21] Fable E., Angioni C., Hobirk J., Pereverzev G., Fietz S. and Hein T. 2011 Confinement and transport properties during current ramps in the ASDEX Upgrade tokamak *Nucl. Fusion* **51** 043006

[22] Fable E. et al 2013 Novel free-boundary equilibrium and transport solver with theory-based models and its validation against ASDEX Upgrade current ramp scenarios *Plasma Phys. Control. Fusion* **55** 124028

[23] Ho A. et al 2023 Predictive JET current ramp-up modelling using QuaLiKiz-neural-network *Nucl. Fusion* **63** 066014

[24] Van De Plassche K.L., Citrin J., Bourdelle C., Camenen Y., Casson F.J., Dagnelie V.I., Felici F., Ho A. and Van Mulders S. 2020 Fast modeling of turbulent transport in fusion plasmas using neural networks *Phys. Plasmas* **27** 022310
[25] Maget P. et al 2022 Healing plasma current ramp-up by nitrogen seeding in the full tungsten environment of WEST Plasma Phys. Control. Fusion **64** 045016

[26] Fietz S., Fable E., Hobirk J., Fischer R., Fuchs C., Pereverzev G. and Ryter F. 2013 Investigation of transport models in ASDEX Upgrade current ramps Nucl. Fusion **53** 053004

[27] Koechl F. et al 2018 Optimising the ITER 15MA DT baseline scenario by exploiting a self-consistent free-boundary core-edge-SOL workflow in IMAS Preprint: IAEA Fusion Energy Conf. (Gandhinagar) p EX/7-25

[28] Coronaa D. et al 2022 Plasma shape control assessment for JT-60SA using the CREATE tools Fusion Eng. Des. **146** 1773-7

[29] Dubbioso S., di Grazia L.E., De Tommasi G., Mattei M., Mele A. and Pironti A. 2019 Vertical stabilization of tokamak plasmas via extremum seeking IFAC J. Syst. Control **21** 100203

[30] de Vries P.C. and Gribov Y. 2019 ITER breakdown and plasma initiation revisited Nucl. Fusion **59** 096043

[31] Imbeaux F. et al 2015 Design and first applications of the ITER integrated modelling & analysis suite Nucl. Fusion **55** 123006

[32] Peeters A.G., Camenen Y., Casson F.J., Hornsby W.A., Snodin A.P., Strintzi D. and Szepesi G. 2009 The nonlinear gyro-kinetic flux tube code GKW Comput. Phys. Commun. **180** 2650-72

[33] Angioni C., Peeters A., Pereverzev G., Ryter F. and Tardini G. (ASDEX Upgrade Team) 2003 Density peaking, anomalous pinch and collisionality in tokamak plasmas Phys. Rev. Lett. **90** 205003

[34] Moret J.M., Duval B.P., Le H.B., Coda S., Felici F. and Reimerdes H. 2015 Tokamak equilibrium reconstruction code LIUQE and its real time implementation Fusion Eng. Des. **91** 1-15

[35] Marin M., Citrin J., Giroud C., Bourdelle C., Camenen Y., Garzotti L., Ho A. and Sertoli M. 2023 Integrated modelling of neon impact on JET H-mode core plasmas Nucl. Fusion **63** 016019

[36] Staebler G.M., Howard N.T., Candy J. and Holland C. 2017 A model of the saturation of coupled electron and ion scale gyrokinetic turbulence Nucl. Fusion **57** 066046

[37] Ho A., Citrin J., Bourdelle C., Camenen Y., Casson F.J., Van De Plassche K.L. and Weisen H. 2021 Neural network surrogate of QuaLiKiz using JET experimental data to populate training space Phys. Plasmas **28** 032305

[38] Romanelli G. and Cenacchi M. JETTO manual (Italy: JET Joint Undertaking) 2.1988, 2

[39] Porcelli F., Boucher D. and Rosenbluth M.N. 1996 Model for the sawtooth period and amplitude Plasma Phys. Control. Fusion **38** 2163-86

[40] Angioni C., Gamot T., Tardini G., Fable E., Luda T., Bonanomi N., Kiefer C.K. and Staebler G.M. 2022 Confinement properties of L-mode plasmas in ASDEX Upgrade and full-radius predictions of the TGLF transport model Nucl. Fusion **62** 066015

[41] Wensing M. et al 2021 SOLPS-ITER validation with TCV L-mode discharges Phys. Plasmas **28** 082508

[42] Bourdelle C., Garbet X., Singh R. and Schmitz L. 2012 New glance at resistive ballooning modes at the edge of tokamak plasmas Plasma Phys. Control. Fusion **54** 115003