Predictive JET current ramp-up modelling using QuaLiKiz-neural-network

This work applies the coupled JINTRAC and QuaLiKiz-neural-network (QLKNN) model on the ohmic current ramp-up phase of a JET D discharge. The chosen scenario exhibits a hollow Te profile attributed to core impurity accumulation, which is observed to worsen with the increasing fuel ion mass from D to T. A dynamic D simulation was validated, evolving j, ne , Te , Ti , n Be, n Ni, and n W for 7.25 s along with self-consistent equilibrium calculations, and was consequently extended to simulate a pure T plasma in a predict-first exercise. The light impurity (Be) accounted for Z eff while the heavy impurities (Ni, W) accounted for P rad. This study reveals the role of transport on the Te hollowing, which originates from the isotope effect on the electron-ion energy exchange affecting Ti . This exercise successfully affirmed isotopic trends from previous H experiments and provided engineering targets used to recreate the D q-profile in T experiments, demonstrating the potential of neural network surrogates for fast routine analysis and discharge design. However, discrepancies were found between the impurity transport behaviour of QuaLiKiz and QLKNN, which lead to notable Te hollowing differences. Further investigation into the turbulent component of heavy impurity transport is recommended.


Introduction
The prediction of temporal plasma evolution is essential for modern fusion devices, in order to estimate and optimize the performance of any given plasma scenario. The demand for this capability is expected to increase as the devices scale in cost, energy output, and complexity, as is the case with ITER. At the start of every tokamak plasma discharge, the plasma is initiated via a breakdown and then a burn-through phase. Afterwards, a set of control trajectories set up the desired steady-state magnetic configuration following the initiation of the plasma current, called the current ramp-up phase. Accurate modelling of this phase is required for detailed pulse design and preparation of ITER scenarios, starting from the prefusion-power-operation (PFPO) campaigns [1]. The lessons learned from the WEST device [2] highlight the necessity and the challenges of successfully developing the current rampup phase, from both an operational and a plasma physics perspective. An examination of the validity of integrated plasma models in this phase helps to identify any gaps in the existing code base and enables these tools to assist in the development and interpretation of these campaigns [3].
Integrated plasma transport models interconnect multiple independent physics models to consistently evaluate the plasma state under different concurrent phenomena. Each of these component models are responsible for computing a specific plasma phenomenon, e.g. transport fluxes, sources or sinks. This study focuses on the JINTRAC [4] model, which contains a time evolution calculation scheme. This allows it to take advantage of the timescale separation between various physics phenomena to optimize and orchestrate the execution of the modules without compromising the physical plasma description.
High-fidelity physics models are often precluded from high throughput analysis workflows due to their heavy computational costs, especially when using a time evolution scheme such as JINTRAC. Within the current tokamak plasma modelling landscape, the simulation bottleneck is typically the turbulent transport calculations. Neural network (NN) surrogate models of the reduced-order gyrokinetic turbulent transport model, QuaLiKiz, aptly named QuaLiKiz-neuralnetwork (QLKNN) [5], provided a promising method of improving this speed without sacrificing as much accuracy as other existing simplified models. The QLKNN family of models combine careful data curation and threshold-aware cost functions to more accurately capture the crucial salient features of plasma microturbulent transport physics. This study focuses primarily on the latest version trained based on JET experimental profile data, labelled QLKNN-jetexp-15D [6], mostly to verify its applicability to a variety of JET discharges.
The coupled JINTRAC plasma transport and QLKNN turbulent transport surrogate model was applied to model a reference deuterium (D) hybrid ramp-up discharge, JET#97776. The hybrid scenario was developed to access a larger β N than the standard H-mode scenarios by preventing the development of large (1, 1) MHD instabilities, or sawtooth crashes, which occur inside the q = 1 surface [7]. While the plasma conditions typically recover between sawteeth crashes, they were also observed to trigger neoclassical tearing modes outside the q = 1 surface which effectively limited the achievable experimental β N [8]. Since the sawteeth instabilities are strongly tied to the q = 1 flux surface, they can be avoided by applying strong safety factor, q, profile tailoring during the ramp-up phase [9][10][11] to keep q > 1 across the entire plasma radius. Once a desirable q profile is established, it is 'frozen' in place by activating the external heating systems, where the resulting high temperatures lower the plasma resistivity and thus limit the current diffusion [12].
Incidentally, since this JET high-β N plasma regime is typically dominated by ion temperature gradient (ITG) driven turbulence, the q-profile tailoring also improves the turbulent transport characteristics during the high-power phase. This is both due to the inner region with low magnetic shear [13] and the suppression of microturbulent instabilities around ρ tor ≃ 0.7 with higher magnetic shear,ŝ [14]. The turbulence in both of these regions are also conducive to electromagnetic stabilization processes [15,16].
However, there are risks that accompany the performance benefits of this hybrid scenario. The hybrid ramp-up at JET exhibits a local decrease of the electron temperature, T e , within the central core, referred to as temperature hollowing. If this effect is strong enough, the core resistivity changes sufficiently to reverse the current density gradient, ∇j, leading to a reversal of the magnetic shear,ŝ. The extrapolation of previous protium (H) experiments indicate that the temperature hollowing effect would be greater in corresponding tritium (T) experiments [17]. This shear reversal can lead to additional MHD instabilities [18], which are observed to be reliable precursors for plasma disruptions [19,20] if they stop rotating in the plasma or lock. For this reason, any detection of their presence generally triggers an early plasma termination within the tokamak control system, making the ramp-up trajectory design more difficult. This study attempts to use the accelerated highfidelity model to identify the mechanism behind the observed temperature hollowing and assist in the development of a successful hybrid ramp-up for the JET T experimental campaign.
This study also serves a secondary purpose other than facilitating a physics analysis of the JET hybrid ramp-up scenario. As the previous QLKNN-jetexp-15D model validation exercises [6] were performed using plasma profiles which were also used in creating the training dataset, this exercise also tests the performance of the NN on data not explicitly included within the training dataset, i.e. whether it can provide accurate estimations for scenarios that were not explicitly used to guide its development.
As the impurity transport is expected to be crucial in modelling this scenario, it is noted that the QLKNN-jetexp-15D model does not directly predict the impurity particle fluxes, Γ imp . Instead, it assumes that they can be acquired by scaling the electron particle flux predictions, Γ e , according to the following equation: where i denotes any generic ion species, including the main fuel ions. As a quasilinear formulation of the turbulent contribution to the impurity fluxes is still under investigation [21], this ad-hoc implementation allows for the imposition of ambipolarity as an interim constraint. Due to the existence of multiple impurity species in the simulation, this constraint does not yield a unique impurity transport configuration on its own, hence the necessity of stipulating the explicit use of equation (1). While the validity and utility of this ad-hoc implementation is discussed in this report, an investigation of the exact turbulent contribution to the impurity particle flux is outside the scope of this study. Section 2 provides the details of the D hybrid discharge, JET#97776, relevant to this study and section 3 describes the development and validation of its reference simulation. Section 4 describes the extrapolation of the reference D simulation to a H and a T plasma, including a discussion comparing the isotopic trends observed in the simulations against known experimental ones. Section 5 then compares the predictive extrapolation simulations against actual T experiments performed with the assistance of said simulations and a deeper examination of the ad-hoc ambipolarity constraint implementation mentioned above. Finally, a summary is provided in section 6 and comments are made on any potential future work.

Reference JET D hybrid discharge
The chosen D reference discharge, JET#97776, was an ohmic test of the proposed ramp-up approach to the JET hybrid scenario in preparation of the JET tritium campaigns in 2021. The reference discharge was performed with I p = 2.26 MA and B T = 3.43 T inside the current flat-top phase. As this exercise focuses on the current ramp-up phase, the actual values of these parameters change over the duration of the analysis window. Figure 1 shows the time trace of the global plasma and pertinent engineering parameters for this discharge.
In addition to clearly exhibiting the hollow temperature profile under investigation, this particular discharge was chosen for this analysis since the auxiliary heating was intentionally omitted from the experiment. For this specific hybrid ramp-up approach, both the NBI and ICRH systems are typically activated at t = 7 s. The lack of auxiliary heating allows for a more accurate depiction of the plasma state at the planned activation time of the auxiliary heating. Additionally, it allowed the current redistribution to continue until the formation of an unstable MHD (1, 1) mode [22], or sawtooth instability, at t ≃ 8 s. The onset time and radial extent of the first sawtooth crash is largely attributed to evolution of the safety factor, q, profile, specifically the q = 1 surface [23], were not activated during this discharge as it was an ohmic test pulse. The fast current ramp allows faster current penetration into the core to flatten the q profile at the cost of a lower q 0 , while the subsequent current drop ensures a higher q 95 right before the heating systems are activated, maintaining a higher magnetic shear,ŝ, to reduce the NTM drive at higher rational q surfaces. This approach leads to the plasma current 'overshoot' feature around t = 5 s [9]. which is itself determined by the current density profile, j. These experimental conditions both simplify the modelling requirements and allow an ideal environment for a simultaneous validation of the current evolution and diffusion model. Table 1 shows the most important details of the JINTRAC simulation using the QLKNN-jetexp-15D model to predict the turbulent transport coefficients. The simulated plasma start time (t = 1.77 s) was chosen to be near the first reconstructed equilibrium within the discharge with a diverted geometry, i.e. after X-point formation, and the end time (t = 9 s) was chosen to be sufficiently past the onset time of the MHD (1, 1) mode to ensure that it was resolved in the simulation. The input data preparation for the initial and boundary conditions of the model was done using the EX2GK tool using Gaussian Process regression fits [24], with the new addition of electron cyclotron emission (ECE) diagnostics for improved electron temperature resolution in the inner core. Due to the importance of the current diffusion in this scenario, the electron temperature [25], T e , boundary condition was deemed a particularly crucial quantity to estimate accurately.

Reference D simulation
The radiative loss terms in the energy balance equations were expected to play a major role in successfully capturing the temperature hollowing effect. This sink term is likely dominated by the heavier impurity species inside the plasma core. Thus, a reasonable estimate of the initial and boundary conditions (IC and BC) for the impurity profiles were required within this simulation in addition to the T e boundary condition. Due to the lack of a transport barrier at the edge of this L-mode plasma and the non-stationary nature of the ramp-up phase, the impurity concentration within the core were computed using the self-consistent evolution equations in the model. It should be noted that a quantitative comparison of the impurity density evolution capability of this configuration of JINTRAC is not well validated.
In order to simultaneously follow Z eff , P rad , and keep the light impurity dilution within acceptable limits, i.e. n Be /n e ≲ 0.03, it was necessary to include three impurity species within the simulation: one light, one mid-range and one heavy impurity species. Due to the specific materials of the JET plasmafacing components, these were selected as beryllium (Be), nickel (Ni) and tungsten (W), respectively. The initial impurity concentrations and profiles were estimated by adjusting them until the fixed-boundary equilibrium calculation module, ESCO [4], converged while using its recommended settings. This resulted in the relative concentrations of 1.0, 0.05, and 0.001 for Be, Ni, and W respectively.
The impurity ICs were originally taken from an interpretive JINTRAC run where the impurity density profiles were identical in shape to the electron density profile but internally scaled to match a flat Z eff = 1.5, i.e. approximately its measured value at the initial plasma time of the simulation. Then, the time-dependent impurity BCs were adjusted until the time evolution of the effective charge, Z eff , and the total radiated power, P rad , agreed reasonably with the Bremsstrahlung optical spectroscopy and bolometer measurements, respectively. These modifications were accompanied by some moderate heuristic changes in the impurity ICs to maintain profile smoothness, resulting in significantly more hollow impurity IC profiles than the original ones. Figure 2 shows the comparison of the simulated time traces to the experimental 0D values for P rad and Z eff . The effective charge measurement after t = 7.5 s implies a pure hydrogenic plasma, within the uncertainty of the measurement. Although Z eff is expected to be low in this phase of the discharge, it was deemed more reasonable to assume a target of Z eff ≃ 1.15 when deciding the impurity densities at those times 5 . The reference simulation was executed using two initial q profiles, each coming from a different Equilibrium FITting (EFIT) [26] equilibrium reconstruction: • one constrained using experimental pressure profiles (EFTPs); • and one constrained using both the EFTPs and polarimetry measurements of the core magnetic field, (EFTF).
This was done to assess the relative impact of the q profile ICs on the time evolution of the simulation. The main difference between the two q profile ICs are primarily seen in q 0 but they converge as the time evolution progresses, as expected after many current diffusion times have elapsed. The effect on the global plasma parameters is negligible. Thus, no further against measured data (gray line) for JET#97776 from t = 1.77 to 9.0 s, using two different initial q profiles from the pressureconstrained equilibrium reconstruction (EFTP-blue) and the polarimetry-constrained equilibrium reconstruction (EFTF-red). The quantities shown are the line-integrated effective charge (top left), Z eff , the total radiated power (top right), P rad . Note that the blue curve lies directly underneath the red curve.
adjustments were made to the impurity density ICs and BCs between the two q profile cases shown in this section. Figure 3 show the time trace comparison of the volumeaveraged electron temperature, ⟨T e ⟩, the volume-averaged electron density, ⟨n e ⟩, the ITER normalized internal inductance [27], l i (3), and the safety factor at the magnetic axis, q 0 , from the reference simulation and their corresponding processed measurement values inside the JET data repository. Good general agreement is achieved with the experimental measurements, with the minor exception of ⟨T e ⟩ being systematically higher. While it is not inherently concerning, the reason for this will be discussed along with the comparison of the 1D T e profiles.
Since the safety factor on the magnetic axis, q 0 , drops below unity in the simulation, an additional validation metric can be performed to ensure that the reference simulation adequately models the dynamic current redistribution. Figure 4 compares the time-dependence of the q profile at the observed inversion radius and the observed sawtooth onset time, given by the vertical dashed line. As shown in the figure, the arrival of the q = 1 surface at ρ tor = 0.2 within the simulation is ∼500 ms earlier than the first experimentally observed sawtooth crash. While further modelling with MHD codes is required to determine the actual sufficiency of the plasma condition to the formation of sawteeth instabilities, this heuristic comparison is taken to indicate that the neoclassical conductivity used within JIN-TRAC, combined with the improved kinetic profiles, provides an adequate description of the current diffusion within these plasma conditions. Although the results are not shown, additional simulations attempted using only the Spitzer conductivity term [28] within JINTRAC resulted in the q 0 dropping to ∼0.5 at t = 3 s. This implies that the neoclassical formulation is required for accurately modelling these plasma regimes. This appears in contrast with recent current diffusion modelling results of the rampup phase in MAST and JET [29]. Future work is recommended in determining whether the discrepancy between this work and the previous modelling exercise came from a difference in plasma regimes or from improved data handling to better capture the T e BC. against measured data (gray line) for JET#97776 from t = 1.77 to 9.0 s, using two different initial q profiles from the pressureconstrained equilibrium reconstruction (EFTP-blue) and the polarimetry-constrained equilibrium reconstruction (EFTF-red). From top to bottom, the quantities shown are the volume-averaged electron temperature, ⟨Te⟩, the volume-averaged electron density, ⟨ne⟩, the line-averaged density,ne, safety factor at the magnetic axis, q 0 , the safety factor near the separatrix, q 95 , and the normalized internal inductance, l i3 . Figures 5 and 6 shows the kinetic profiles for the electron temperature, T e , and electron density, n e profiles, respectively, at various time slices within the simulation, compared against experimental measurements averaged within a 250 ms window around the stated simulation time. The electron temperature measurements were taken from two separate Thomson scattering diagnostics (LIDR and HRTS) and the ECE diagnostic. The density measurements were taken from the same two Thomson scattering diagnostics as the electron temperature. These figures indicate that the density peaking and temperature hollowing effects observed in the measurements are adequately captured, along with their time evolution. A more detailed discussion of the heat and particle fluxes predicted by the NCLASS neoclassical transport model [30] and the QLKNN turbulent transport model for this simulation can be found in appendix A.    The initial q profile switch altered the depth of the T e hollowing as the plasma evolved, but it is uncertain which q profile IC is more accurate due to the discrepancy at ρ tor = 0.3. Since the initial q profile has only a minor global impact on the simulation and the ESCO equilibrium boundary is taken from the EFTF calculation, the EFTF q profile run was used as the reference simulation for this study. Additionally, it is noted that the simulated q profile does not agree well with the estimate from the EFIT algorithm. This further underlines the necessity to perform a predictive current simulation to selfconsistently relax the current profile before using it within an integrated modelling exercise, as is current recommended practice within the field. Also, the discrepancy between the simulated and measuredT e shown in figure 3 can be attributed to the increased T e around ρ tor = 0.3, although it is unclear if this is a result of the neoclassical transport or the NN turbulent transport predictions at this moment. Figure 7 shows the time traces for the predicted volumeaveraged impurity density concentrations for the three impurity species (Be, Ni, and W) and figure 8 shows the density profiles for each of the three impurity species at the same time slices as above. It can be directly seen from these profiles that the temperature hollowing comes from the accumulation of high-Z impurity species in the central core, where the peaking is much stronger for heavier impurities than light ones. This is the expected behaviour from the neoclassical inward flux for heavy impurities in this plasma scenario [31], although it is worth noting that the turbulent transport contribution to the impurity fluxes do not seem to heavily modify this response.

Role of impurities on observed temperature hollowing
Since the central temperature hollowing effect is tied to the impurity transport, two possible explanations were provided regarding the mechanism: 1. through the increased local effective charge, Z eff , which lowers the neoclassical conductivity [25], preventing current penetration and reducing the ohmic heating term in the energy conservation equation; 2. and through the increased local radiated power density, Q rad , which acts as a sink term in the energy conservation equation.
Note that since these two channels stem from different physics phenomena affecting the same plasma quantity, it is suspected that both are present but it is crucial to understand their relative contribution. Due to the success of the reference simulation on capturing the T e hollowing, additional simulations were performed to attempt to determine the relative importance of these two mechanisms. The main conclusion of this separation study is that Z eff only affects the central T e as an integrated effect from the plasma edge, whereas Q rad has a more local impact. This confirms that an increased T e hollowing within the scenario reflects an increase of heavy impurity accumulation specifically within ρ tor < 0.2, as opposed to an increase only inside the turbulence region of ρ tor ⩾ 0.4. It also implies that the observed temperature hollowing is more dependent on Q rad than the neoclassical conductivity. For brevity, the specific details of this separation study are both shown and discussed further in appendix B.
Additionally, as indicated by table 1, the reference simulation was executed assuming a fixed zero toroidal angular frequency, Ω tor , profile and without including the neutral source code, FRANTIC [32], for the edge particle source. The former was justified as the low levels of toroidal momentum due to the lack of NBI heating in JET L-modes typically have a negligible impact on the turbulent transport calculation. The latter was justified by the relatively inward position of the simulation BC, ρ tor = 0.9, which is beyond the expected region where the ionizing neutral source term plays an important role in the main ion transport equations. However, in the interest of completeness regarding the role of impurities on the T e hollowing effect, two additional independent simulations were performed as sensitivity studies: • one including a non-zero Ω tor profile taken from another ohmic L-mode flat-top discharge, as no rotation measurements were available for JET #97776; • and one including the FRANTIC model with timedependent inputs adjusted to provide the experimental gas puff injection rate.
Both of these simulations resulted in differences of <3% for the predicted n e , T e , and T i profiles compared to the reference simulation presented earlier. As these effects are negligible to the analysis of the T e hollowing effect, these contributions were both omitted from the reference simulation in favour of simplicity and their results not explicitly shown in this paper.
However, it is noted that their impacts on the impurity profile evolution are slightly larger. The core n Ni and n W decreased by ∼10% when including the non-zero Ω tor , likely resulting from the added contribution of the toroidal velocity component of the radial electric field on the neoclassical transport calculation [33]. Also, the overall n Be profile decreased by ≲5% when including the FRANTIC neutral source model, accompanied by a ∼5% drop in the line-averaged Z eff and a similar decrease in core n Ni and n W by ∼10%. While these differences in impurity profiles have a negligible impact on the other kinetic profiles within the analysis of this discharge, this may not apply in a general sense. Thus, it is still recommended that these effects be investigated further when analyzing the impurity transport in the plasma ramp-up regime.

Predict-first T simulations
Since the current ramp-up approach for this scenario has already been developed and tested in the device using the D reference discharge, this predictive study assumes it can serve as an initial guess for the engineering waveforms required in the T discharge. The integrated modelling framework can then be used to investigate the performance of the scenario in tritium and determine which factors are the most effective for configuring the associated paired experiment. Specifically, the T e and q profiles can then be monitored to ensure that the T ramp-up approach follows the reference D discharge as closely as possible to avoid triggering the disruption mitigation system.
With the reference simulation configured and its results verified, as discussed in section 3, incremental modifications can be made to it in order to predict the outcome of future experimental scenarios. Over the course of this study, these modifications include: • switching the main ion isotope; • adjusting the electron density, n e , boundary condition; • and adjusting the impurity density, n Be , n Ni , n W , boundary condition.
For brevity, only the modelling results of the first two items are described in the main body of this paper based on their relevance to the main conclusions. The last item is detailed further in appendix C.

Extrapolation to T plasma
This section examines the effects of switching the fuel mix in the experiment from 100% D to 100% T without adjusting any other simulation input. For a more complete picture of the isotope effect, this exercise was also performed using a 100% protium (H) fuel mix. Figure 9 shows the comparison of the time traces of the global plasma parameters using the different main ion isotopes in the simulation. In general, there is no significant difference in the global plasma parameters with respect to the isotopic mass. Some minor effects are the increase of total P rad with isotope mass after t = 7 s without a corresponding increase in global Z eff , indicating a higher heavy impurity content with isotope mass, as discussed in section 3.1. Figure 10 shows the kinetic profiles only at t = 7 s, as it represents the plasma state when auxiliary heating is switched on to effectively freeze the q profile. Figure 11 shows the impurity density profiles at this same snapshot, t = 7 s. From these figures, the main differences between the different simulations indicate that an increase in the main fuel ion mass leads to: • an increase in T e hollowing; • a decrease in T i across the entire radial profile; • an increase in the n e peaking; • a delay in the arrival of q = 1 surface at ρ tor = 0.2; • an increase in the core accumulation of heavy impurities, Ni and W; • and an increase of n Ni and n W around ρ tor ≃ 0.5.  While the T e and q impacts are expected trends from previous experiments comparing hybrid ramp-up phases in H and D plasmas [17], the T i and n e are not. The T i trends may have been present in the previous isotope experiments but the measurements necessary to determine this were not available due to the lack of NBI injection in these discharges.
The lower T i within the T plasma simulation is likely an effect of the reduced electron-ion energy exchange due to the increased main ion mass. It is postulated that this lower T i is the source of the increased temperature hollowing, as it reduces the neoclassical temperature screening effect. Under these conditions, this translates into an increase in the inward flux of impurities leading to increased core impurity accumulation. This then increases the radiated power density, Q rad , in the central core which is responsible for the T e hollowing in that region, as detailed within appendix B. This lower core T e drives the collisionality up, leading to a local increase in current diffusion and generating the reverse shear in the core region if the T e hollowing is severe enough.
To further emphasize the dominant impact of the local radiated power density, Q rad , on the temperature hollowing over the neoclassical conductivity, figure 12 compares the local ohmic heat source, Q ohm , the local radiation power density, Q rad , and the local electron-ion energy exchange source, Q ex across the isotope scan. As mentioned in section 3.1, one expected result of the increased core collisionality is the reduction of current diffusion via a lower local neoclassical conductivity, resulting in a lower central Q ohm . However, figure 12 shows that the central Q ohm actually increases with isotopic mass. Since the ohmic heating process itself occurs via collisional processes, it is suspected that the lower current diffusion is more than compensated by the increased collisional energy transfer. In this case, the temperature hollowing truly only occurs due to the significant increase in Q rad incurred by the increased influx of heavy impurities.
Regarding the other simulation results, the increased delay of the q = 1 surface arrival with the isotope mass is an experimentally observed phenomenon [17], but the time difference is quantitatively larger within the experiments. Also, the increasing n e peaking prediction within the simulation goes contrary to the experimental trends, but the quantitative magnitude of this trend is sufficiently small to be considered negligible. Due to the contribution of n Be to n e , it is also possible that the reversal of the trend is due to quantitative discrepancies in the impurity transport coefficients between the model and experiment. A few additional insights into these discrepancies are provided in section 5.

Advising the development of the paired T experiment
Since both current diffusion and ohmic heating are modified by the plasma collisionality, it is suspected that the q and T e profiles can be affected by changing the plasma density, indicated by n e . This fact has been experimentally exploited in previous H discharges to construct a paired experiment with the same q profile as the D reference discharge. Comparing these experiments allow the separation of phenomena which are affected directly by the main ion mass from those which are affected indirectly via another plasma parameter.
Within the previous JET H discharges, the core plasma density was controlled by injecting neutral gas into the SOL. While other methods are available for controlling the core impurity accumulation [2,34,35], this was ultimately chosen due to experimental constraints but remains suitable for this study due to its simplicity in modelling. By empirically extrapolating the gas flow rate necessary to produce the paired experiments in H [17], it was estimated that an increase of 20%-30% was needed in the gas injection system in order to achieve the same q profile in T. This range accounts for the empirical uncertainties regarding the effect of impurities and other isotope effects. This modelling study was performed to justify and/or refine that empirical extrapolation, as well as provide key physics insights into the impurity and other isotope effects which affect that estimate.
As mentioned earlier, the integrated model settings chosen for this exercise do not include any regions beyond the internal boundary condition, including the separatrix region and scrape-off layer (SOL). Thus, the effect of the additional gas injection cannot be directly modelled but a proxy can be made by adjusting the time-dependent n e BC (ρ tor = 0.9) in the simulation inputs. Figure 13 shows the results of the n e BC scan on the kinetic profiles at t = 7 s,  including the self-consistent evolution of the density profile. Figure 14 shows the results of the n e BC scan on the impurity density profiles at t = 7 s and figure 15 shows the results of the n e BC scan on the major heat source and sink profiles.
As expected from the H experiments, the higher n e from the BC modification results in less T e hollowing. This is due to the increase in the plasma collisionality, which adjusts the turbulent transport such that it lowers the overall density gradient in the profile. This in turn reduces the neoclassical inward flux of heavy impurities, leading to less impurity accumulation and a lower Q rad in the central core region. Due to the explicit connection of this mechanism to the ITG turbulence regime present in this scenario, this method of controlling T e hollowing via gas puff is likely not universally applicable.

Comparison with post-analysis T experiment
A pair of discharges were performed incorporating the results of the predictive modelling exercise discussed above. This pair of discharges are hybrid ramp-up discharges with a similar approach to the reference discharge (JET#97776) but using tritium as the main fuel. One (JET#98562) was performed at the reference gas flow rate, which ultimately led to a triggering of the disruption mitigation system, and the other (JET#98567) was performed with an increased gas flow rate targeting a ∼20% higher line-integrated density, which did not trigger the disruption mitigation system. It is noted that the exact percentage increase varied with time in the T experiment although this set a valuable target.
While the success of the predictive modelling exercise in quantitatively advising target control room parameters is already a milestone achievement for the modelling suite, it is also interesting to inspect the quality of the predicted 1D profiles. Figure 16 shows the comparison of the predicted kinetic profile from extrapolating the D reference to T against the experimental measurements taken from this pair of tritium discharges.
From these figures, the n e and T e predictions between ρ tor ≃ 0.3 and the simulation boundary shows excellent agreement with the experimental measurements within their uncertainties. However, discrepancies beyond those uncertainties are found in the central region, from 0 ⩽ ρ tor ≲ 0.3. This implies that the impurity profile predictions, and consequently the impurity transport coefficients from QLKNN, may not be fully representative of the physical scenario. A deeper investigation into the turbulent component of the impurity transport Figure 16. Comparison of the predictive T extrapolation simulations from the D reference discharge (solid lines) and recent JET T experiments (gray points). The scenario using the same absolute gas flow as the D reference (JET#98562-top) resulted in a mitigated disruption, whereas the scenario with an increased gas flow (JET#98567-bottom) did not. The profiles agree for ρtor ≳ 0.3, but significant discrepancies occur within the central core.
behaviour is recommended, but remains outside the scope of this study.

Comparison with QuaLiKiz model
In order to complete the QLKNN validation exercise, it is prudent to compare the performance of the reference discharge using QLKNN against the original QuaLiKiz model. However, this comparison is complicated by to the impurity transport assumption made in implementing the QLKNN within JINTRAC.
Firstly, an initial comparison was made including a modification to the QuaLiKiz implementation within JINTRAC such that it also used the ad-hoc impurity transport scaling described by equation (1). Figures 17 and 18 show the comparison between these two simulations, where all the QLKNN kinetic profiles fall within the ±10% profile RMS error with respect to QuaLiKiz, as quoted by [6].
These plots confirm that the QuaLiKiz electron and main ion transport coefficients in this physical scenario are accurately replicated by the QLKNN-jetexp-15D model, effectively providing additional validation for those output variables. With this confirmed, the next step is removing the ad-hoc impurity transport modification to determine whether it is the source of the discrepancy and/or the extent of its impact on the simulation as a whole. Figure 19 shows the results of the D Figure 17. Profile comparisons of the electron temperature (left), Te, ion temperature (center), T i , and electron density (right), ne, using QLKNN (blue line) and a modified QuaLiKiz implementation (red line) as the turbulent transport model. The QuaLiKiz implementation in JINTRAC was modified by replacing the computed impurity particle flux, Γ imp , with a scaled electron particle flux, Γe, according to equation (1). Figure 18. Profile comparisons of the Be impurity density (left), n Be , Ni impurity density (center), n Ni , and W impurity density (right), n W , using QLKNN (blue line) and a modified QuaLiKiz implementation (red line) as the turbulent transport model. The QuaLiKiz implementation in JINTRAC was modified by replacing the computed impurity particle flux, Γ imp , with a scaled electron particle flux, Γe, according to equation (1). This shows that the original QuaLiKiz impurity profile predictions vary significantly from those of the ad-hoc rule. The increased impurity density in the central core, specifically the radiating species of Ni and W, explains the increased T e hollowing. Additionally, the variations in the impurity density and its gradient may alter the turbulent transport characteristics sufficiently to change the n e profile, especially near the simulation boundary.
A more direct and quantitative comparison of the Ni and W particle flux over the initial time steps of the simulation is provided in figures 20 and 21 respectively. This indicates that the ad-hoc ambipolarity rule assigns an outward flux across the predictive region in the first 100 ms, which helps to establish a degree of hollowness in the impurity density profile. Then, In contrast, the QuaLiKiz model provides a much higher inward flux from the beginning of the discharge and remains heavily inward as the discharge progresses. This leads to a notably higher density at ρ tor = 0.2, which is then transported into the central core via neoclassical transport. This leads to a much higher impurity accumulation and a more pronounced T e hollowing effect. While the turbulent component of impurity transport has not been fully validated within the QuaLiKiz code, the importance of its accurate prediction in the ramp-up phase is highlighted within this study. While a deeper fundamental investigation into this phenomenon is outside the scope of this work, this paper strongly recommends it before these reduced models are used for extensive prediction in non-Hmode plasma regimes.

Conclusion
Overall, the JINTRAC integrated model, coupled with the newly-developed QLKNN-jetexp-15D model for turbulent transport predictions, was successfully validated on the deuterium hybrid ramp-up discharge, JET#97776, and further used for predictive simulations for tritium plasma operation at JET. This predictive capability of the model is verified not only via the simulation of the plasma ramp-up regime, a novel application for the QuaLiKiz model in general, but also the via the extrapolation of the scenario to a pure T discharge, returning results which are in good agreement with isotopedependent trends from previous H discharges. With T, the T e hollowing increases and the onset of the (1, 1) MHD sawtooth instability may be delayed as indicated through the q = 1 arrival time. The only trend which appears in opposition to the experimentally observed trends is the density peaking, which increases with isotope mass in the simulations but decreases in experiments. However, the relative change of the density peaking from D to T is small enough in both simulation and experiment that both fall within the GPR fit uncertainties and may be considered negligible.
Aside from the demonstration of the predictive capabilities of the integrated modelling suite, this study also provides some insight into the physical process causing the increased temperature hollowing observed when moving from H to D in the previous isotopic experiments. The simulations in this study indicate that the increased main ion mass lowers T i through a lower electron-ion heat exchange. This increases the neoclassical inward flux of the heavy impurities (e.g. Ni and W in this study) and leads to more accumulation in the central core. These heavy impurities then act as a local heat sink via radiation losses, effectively lowering the core T e . The consequence is then a higher collisionality, increased current diffusion, and a higher q 0 . While each of these individual mechanisms and correlations are well-documented within the literature of their respective specialized physics domains, the chain of interactions only became apparent through the analysis of this scenario within an integrated model. A key advantage of this approach is the capability to determine the relative strength of a given interaction compared to all the others present in a given scenario, as evidenced by this study.
The validation of the QLKNN-jetexp-15D simulations against the original QuaLiKiz model was complicated by the use of the ad-hoc ambipolarity constraint for the impurity flux, as QLKNN-jetexp-15D was not trained directly to predict turbulent impurity transport coefficients. In terms of the electron and main ion coefficients, the two models show good agreement as shown using an implementation of QuaLiKiz which uses the same ad-hoc ambipolarity constraint. However, without this constraint, the two models exhibit significant discrepancies which are attributed to the differences in the turbulent predictions of the impurity transport coefficients. A deeper fundamental study into the validity of these impurity transport coefficients within the original Qua-LiKiz model are recommended, since the predicted profiles with the original model also deviated from the experimental measurements.
Finally, the increased computational speed of the QLKNNjetexp-15D model was leveraged to assist with answering operational questions within the time constraints of the planned experimental schedule. Each full ramp-up simulation took ∼6 h on a single CPU using the NN, whereas the original QuaLiKiz code took ∼130 h on 16 CPUs. This allows not only for the rapid development of a state-of-the-art validated rampup simulation reference but also for the execution of numerous sensitivity studies required to provide meaningful physics insights. Approximately 500 simulations were completed over the span of 4 calendar weeks for the bulk of this study. In the end, the predictive study estimated a required gas injection rate increase of 20%-30%, corresponding to a line-averaged density increase of ∼20%, with reasonable assumptions on the impact of the gas injection on the plasma edge and SOL. This agrees the empirical trends from isotope effect studies from previous H discharges, which was extrapolated to T for an estimated gas injection increase of 30%. The effectiveness of this qualitative recipe was further corroborated by a pair of JET T discharges and provides even greater confidence in the predictive capabilities of the JINTRAC integrated model, using the QLKNN-jetexp-15D turbulent transport predictions, within this plasma regime at JET.

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) and from the EPSRC (Grant Number EP/W006839/1). 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. Neither the European Union nor the European Commission can be held responsible for them. Figures 22 and 23 show the heat fluxes, q, and particle fluxes, Γ, respectively, for various species in the plasma at specific time slices of the reference simulation. It is important to restate that the turbulent contribution to the ion particle fluxes is proportional to the electron particle flux, following the ad-hoc relation described by equation (1). Thus, all the profiles of the turbulent contributions in the reference simulation, shown in figure 23, are the same and the only variation between the species is the neoclassical contribution. A rudimentary estimation of the particle confinement time, τ p , indicated by this level of turbulent transport can be computed using the following expression:

Appendix A. Reference simulation transport coefficients
where n s is the number density of species s, V is the volume contained within the flux surface located at the given ρ coordinate. This calculation yields a particle confinement time on the order of ∼1 s for the QLKNN simulation, which is consistent with neoclassical transport time scales. Deeper gyrokinetic studies on turbulent impurity transport in this regime are required to determine whether the transfer of electron particle transport coefficients, which yields a reasonable simulation of the T e hollowing, more accurately describes the physical scenario or was a fortuitous result.  . As this ramp-up plasma is not in steady-state (dn/dt), it is not expected that these plots show zero total particle flux despite having no internal particle sources.
Concerning the discrepancies observed in section 5.2, it appears that the turbulent component of the impurity particle flux in the original QuaLiKiz model can dominate over the neoclassical component, whereas this appears to not be the case when using QLKNN-jetexp-15D. Note that the comparison is made here on the particle flux, Γ, which already combines the contributions based on their diffusive and convective coefficients, D and V respectively, according to the following expression: for a given species, s, and its number density, n s . This does not imply that the neoclassical diffusive and convective coefficients are of the same order of magnitude as their turbulent counterparts. These separated diffusive and convective components are not shown, as the particular scheme chosen to transfer QLKNN predicted quantities to the impurities, provided in equation (1), renders the separated coefficients meaningless without additional assumptions. When comparing QLKNN-jetexp-15D to QuaLiKiz, a sign change is observed in the Ni and W impurity particle fluxes which leads to significantly different T e hollowing dynamics in their respective ramp-up simulations. One possible explanation is the inadequacy of the electrostatic assumption of Qua-LiKiz to model turbulent impurity transport, as exampled by studies in high-β regimes [36].

Appendix B. Separation of effective charge and radiation effects
In order to separate the effects of Z eff and Q rad on the temperature profiles, the time evolution of the n e and all n imp were imported from the reference simulation and fixed within a new one. This allowed an independent variation of the impurity density profiles for this investigation, where: • Z eff was isolated by modifying the Be impurity profile, since its Q rad contribution is negligible at these core plasma temperatures; • and Q rad was isolated by modifying the W impurity profile, since its Z eff contribution is minor due to its low absolute density in the plasma.
By investigating these effects through the adjustment of the impurity densities, the self-consistent calculation capabilities of the integrated modelling framework can be used to develop a more physically tangible intuition within the study. In this case, this is done via the SANCO [4] impurity transport module. Although the impurity density evolution equations are not solved in these simulations, the module is still used to compute the radiated power based on the charge state distribution and the ADAS cross-section library. Figure 24 shows the modified impurity density profiles used in these simulations along with their effects on the Z eff and P rad time traces. The modifications were added in two ways: Figure 24. Example impurity density profile modifications applied to the prescribed integrated model inputs (top row) for separating local effective charge, Z eff , effects (via Be-left column) and local radiated power density, Q rad , effects (via W-right column). The impact of these modifications on the simulated time traces of total radiated power, P rad , (middle row) and the line-integrated effective charge, Z eff , (bottom row) are also provided for demonstrate its validity.
• using a constant multiplication factor across the entire profile and a modification applied only in the inner core; • and using a Gaussian multiplication factor centred on ρ tor = 0 with σ = 0.067.
These two different modification methods help to distinguish which effects are localized to the inner core region and which are effectively integrated across the entire plasma radius. Figures 25 and 26 show the sensitivities due to changes in the local Z eff , via Be profile modifications, and the local Q rad , via W profile modifications, respectively. Only the profiles at t = 7 s are provided since the effect over the time evolution is fairly consistent.
From these figures, the temperature hollowing is affected more by the local Q rad channel than the local Z eff channel in the central core region. Also, the Q rad contribution can be considered a local effect as the uniform multiplication factor has the same impact on the inner core temperature as the Gaussian multiplication factor. On the other hand, the Z eff contribution is an accumulated global effect as the uniform multiplication factor shows a continuous divergence of the profile from the Figure 25. Comparison of the electron temperature (left), Te, ion temperature (center), T i , and safety factor (right), q, profiles for the local effective charge, Z eff , sensitivity study of JET#97776 at t = 7 s. The dashed vertical lines indicates the internal simulation boundary, ρtor = 0.9, outside which the profiles are prescribed in the simulation. The impact of the Gaussian modification factor is almost negligible but the effect integrated across the whole plasma is noticeable when the uniform multiplication factor is applied. Figure 26. Comparison of the electron temperature (left), Te, ion temperature (center), T i , and safety factor (right), q, profiles for the local radiated power density, Q rad , sensitivity study of JET#97776 at t = 7 s. The dashed vertical lines indicates the internal simulation boundary, ρtor = 0.9, outside which the profiles are prescribed in the simulation. The impact across the entire radial profile is almost negligible but the effect on the central core is roughly equivalent regardless of the uniform multiplication factor or the Gaussian multiplication factor. edge while the Gaussian multiplication factor has nearly no impact.
This Z eff modification result also raises the possibility that the increased local effective charge, Z eff , also lowers the turbulent transport coefficients within the central core. This can be done either through a stabilization of TEM modes by increasing the collisionality, or through a stabilization of ITG modes due to the impact of light impurity density and/or density peaking [37]. This is further supported by the fact that the q profiles for all the various Z eff variations are nearly identical, although the plots for this are not provided in this report for brevity.

Appendix C. Impurity boundary condition scans
Due to the higher mass of tritium over deuterium, an increase in the sputtering yield from the divertor and other plasma facing components is expected with the switch of the main ion [38], assuming all other factors remain the same. Since these simulations do not extend into the plasma edge or SOL regions, a quantitative estimate of the sputtering source and the impurity transport into the plasma core was not selfconsistently calculated. As a proxy, the impurity density boundary conditions was varied to study the sensitivity of the core plasma conditions to the increased sputtering yield. Figure 27 show the effect of an impurity density BC scan on the kinetic profiles of the tritium simulation at t = 7 s. The BC scan increases all the impurity BCs, i.e. for Be, Ni, and W simultaneously, by a factor, ∆n imp,m /n imp , of 1.2 and 1.5. The time traces are not shown as they do not add much value to this discussion.
As expected from the fact that the plasma profiles are not significantly altered, the amount of core impurity accumulation increases as the impurity BC increases, leading to more T e hollowing. However, as both T e and T i are also slightly affected for ρ tor ≲ 0.5, this core increase is not necessarily linear with the BC. This effect on the temperature profiles between 0.2 ⩽ ρ tor ≲ 0.5 is likely due to the increasing Z eff with higher impurity content. This can increase the critical threshold of the turbulent transport, depending on the unstable modes present in the scenario, reducing the turbulent transport and allowing the energy conservation equation to support higher gradients. In spite of this, the boundary impurity density interestingly does not seem to have a strong impact on n e , only minimially increasing it within the inner core. Figure 27. Comparison of the electron temperature (leftmost), Te, ion temperature (left middle), T i , electron density (right middle), ne, and safety factor (rightmost), q, profiles at t = 7 s for the impurity BC scan of the T extrapolation simulation (blue) based on JET#97776, using a multiplication factor of 1.2 (red) and 1.5 (purple). The dashed vertical lines indicates the internal simulation boundary, ρtor = 0.9, outside which the profiles are prescribed in the simulation. All impurity BCs were multiplied by the same factor and was done to represent the increased sputtering yield from the SOL due to an increased main ion mass.
Similarly to the switch from 100% D to 100% T, the gas injection rate used to modify the core density also has an impact on the impurity sputtering. This is especially true within the divertor region [39]. The increased neutral density causes the plasma to cool faster from collisional and radiative processes on its approach to the divertor plates. This lowers the incident energy of the particles striking the solid surface and lowers the impurity sputtering rate, specifically of W in the case of the JET divertor plates.
In order to capture this, the time-dependent density profiles were again prescribed and a multiplication factor was applied to all of the impurity density BCs simultaneously. As the exact relation between the gas injection rate and the sputtering yield is not quantitatively known, this study assumes a proportional relation for the impurity density BC modifications in order to this study the sensitivity of this process, given as follows: where c Y = 0.5. For example, this means a 50% increase in the n e BC results in a 25% decrease in the n imp BCs. While it is unlikely that c Y is a constant across all density modifications and all impurity species in reality, this assumption is sufficient for this sensitivity exercise. Figure 28 shows the results of the n e BC scan on the base simulation using prescribed density profiles and figure  29 shows the results of the n e scan on the major heat source and sink profiles. Surprisingly, using the plasma conditions of this scenario, the effect of the combined electron and impurity density BC scan is nearly imperceptible from the results shown in figure 13, except for a marginally lower T i .

Appendix D. Estimated gas injection modification
In order to use these simulations to provide a prediction for the required gas injection, the BC modifications performed to generate figures 13 and 28 can be permuted to generate a grid of simulation results. Then, since the performance of the highpower phase of the hybrid scenario is dependent on the q profile when the auxiliary heating systems are switched on, its properties at t = 7 s can be used to compare the T simulation configuration to the reference D simulation. Specifically, the safety factor at the magnetic axis, q 0 , at t = 7 s is used as the metric since it is the location with the greatest variability in the set of simulated q profiles. Table 2 contains the values of q 0 at t = 7 s for the various permutations of the density BCs used in this predictive study. The target value, q 0 = 1.132, was taken from the validated reference D simulation. It is suspected that the impurity scenarios with ∆ n imp,m /n imp ⩾ 20% and c Y = 0.5 are the most realistic, resulting in the required density increase at the simulation BC to be estimated at 30%-40%. However, the experimental deductions were all made relative to a line-integrated density measurement, a more meaningful comparison would be to compute the line-integrated density of the simulation via a synthetic diagnostic. This is because the self-consistent calculation of the density profile evolution within the original n e BC scan can affect the normalized density gradients such Figure 29. Comparison of the local ohmic heat source (left), Q ohm , local radiated heat source (center), Q rad , and local electron-ion energy exchange source (right), Qex, at t = 7 s for the combined electron and impurity density BC scan of the T extrapolation simulation (blue) based on JET#97776. The plotted region ends at the internal simulation boundary, ρtor = 0.9, due to the large non-physical values present beyond this region.  that the whole profile is not just a scaling of the reference simulation. Figure 30 shows the results of the synthetic diagnostic on the simulation results, as applied to the last three columns of table 2. From these plots and the estimated range of n e BC increase given above, the synthetic diagnostic estimates the corresponding line-integrated n e increase to be between 20%-30%. This value is slightly lower but agrees well with the empirical estimation from the previous H ramp-up experiments, providing an extra degree of confidence in the predictive capability of the JINTRAC and QLKNN-jetexp-15D coupled model within JET ramp-up scenarios.

Appendix E. Inclusion of flux compression via dB 0 /dt
As indicated in figure 1, the toroidal magnetic field, B 0 , is still varying in time for the initial 5 s of the analysis window. This adds an additional complication to the JINTRAC modelling, as the resulting flux compression can both alter the time evolution of the pressure profile and the magnetic geometry. While this impact was assumed to be minor for the conclusions derived in this study, it was deemed prudent to generate a custom build of JINTRAC was developed to include the dB 0 /dt-terms in the current diffusion equation to verify this assumption.
As shown in figure 31, this flux compression term effectively increases the current diffusion in the earlier phases of the simulation, leading to a lower q 0 at t = 7 s and a faster approach of the q = 1 to the inversion radius-about 500 ms faster than the reference D simulation. The compression of the flux surfaces increases the core pressure, via both the density and temperature, while also pushing impurities inward and accelerating the core impurity accumulation, but an overall ∼10% higher core T e and ∼5% higher core n e at t = 7 s. However, this results from a ∼10% increase in Q ohm , a ∼25% increase in |Q rad |, as shown in figure 32. In this scenario, the combined impact of the increased current diffusion and the core pressure increase from the flux surface compression more than  Comparison of the local ohmic heat source (left), Q ohm , local radiated heat source (center), Q rad , and local electron-ion energy exchange source (right), Qex, at t = 7 s for the combined electron and impurity density BC scan of the T extrapolation simulation (blue) based on JET#97776. The plotted region ends at the internal simulation boundary, ρtor = 0.9, due to the large non-physical values present beyond this region.
compensates for the higher radiative loss due to the likewise compression of the impurity species densities.
As the impact on the predicted profiles is mostly limited to the central core region, the simulated profile sensitivities from the dB 0 /dt-terms do not have any strong impact on the conclusions made by this study. However, its contribution to the overall discrepancy in the central core is not negligible and should be considered an important inclusion in future ramp-up scenario modelling. That said, the likelihood for power reactors to employ a ramp-up scenario containing a time-dependent B 0 remains uncertain.