Scenario optimization for the tokamak ramp-down phase in RAPTOR: Part B. safe termination of DEMO plasmas

An optimized plasma current ramp-down strategy is critical for safe and fast termination of plasma discharges in a tokamak demonstration fusion reactor (DEMO), both in planned and emergency scenarios, avoiding plasma disruptions and excessive heat loads to the first wall. Plasma stability limits and machine-specific technical requirements constrain the stable envelope through which the plasma must be navigated. Large amounts of auxiliary heating are required throughout the ramp-down phase, to avoid a radiative collapse in the presence of intrinsic tungsten and seeded xenon impurities, as quantitatively estimated in this work. As the plasma current is reduced, the current density becomes increasingly peaked, reflected by a growing value of the internal inductance ℓi3 , resulting in reduced controllability of the vertical position of the plasma. The feasibility of different plasma current ramp-down rates is tested by applying an automated optimization framework embedding the RAPTOR core transport solver. Optimal time traces for plasma current Ip(t) and plasma elongation κ(t) are proposed, to satisfy an Ip -dependent upper limit on the plasma internal inductance, as obtained from vertical stability studies using the CREATE-NL code, as well as a constraint on the time evolution of q 95, to avoid an ideal MHD mode. A negative current density near the plasma edge is observed in our simulations, even for the most conservative Ip ramp-down rate, indicating significant transient dynamics due to a large resistive time.


Introduction
Safe termination of a burning plasma in a DEMO reactor is a highly non-trivial task.The high energy content of a reactorgrade plasma and the limited thickness of a reactor first wall (to allow for sufficient tritium breeding), lead to a very small tolerance for plasma disruption events.As safe termination scenarios, both for routine and for emergency shutdown, are critical to any viable tokamak reactor concept, present-day tokamak experiments have started to investigate stable ramp-down solutions, guided by modeling tools.
Modeling and optimization with the RAPTOR code, applied to TCV, ASDEX Upgrade (AUG) and JET, has allowed to increase the plasma current ramp-down rate while maintaining the plasma radially and vertically stable, by optimizing the time evolution of plasma current and elongation throughout the ramp-down phase [1].On DIII-D, vertical displacement events (VDEs) could be successfully prevented by adjusting the plasma elongation and the inner-gap to the vessel wall, in response to real-time estimators of the proximity to the vertical stability limit [2].Furthermore, [2] presents the analysis of a large data set of emergency shutdowns after large tearing modes.Modeling of the ramp-down phase for AUG is discussed in [3,4] and for JET in [5][6][7].
Fast simulators like RAPTOR allow to systematically explore reactor operating points and scenario dynamics, and to automatize the search for optimal control strategies.In [8], METIS [9] is used to explore operating points for a pulsed and a steady state DEMO design.Note that the large size and the high temperature of a DEMO plasma core lead to a very slow diffusion time scale of the current density.This clearly illustrates the need for a fast simulation tool to optimize a DEMO scenario.
While we focus on ramp-down optimization for a DEMO plasma, similar considerations are relevant for ITER.Timedependent, self-consistent ramp-down simulations are needed to project how observations on present-day devices will scale to ITER operation, due to non-linear dynamics and the varying characteristic time scales at play, as argued in [10].The complexity of this task, requiring simultaneous accounting for plasma position and shape control, energy and density evolution, MHD mode activity and impurity transport, demands an analysis that combines a range of time-dependent equilibrium and transport solvers of varying fidelity.The increase of internal inductance after the HL transition and the need for a significant elongation reduction during ramp-down, resulting in an ITER termination scenario remaining at relatively low q 95 ∼ 3 have been discussed in [11].Furthermore, [11] makes the observation that I p will be an effective control parameter for the internal inductance evolution, as the ITER ramp-down is fast with respect to the resistive time τ R .Strategies for pellet fueling and auxiliary heating to avoid tungsten accumulation during the ITER baseline ramp-down have been studied in [7,12].
In the previous paper, 'Part A: Analysis and model validation on ASDEX Upgrade' [13], we have discussed how RAPTOR has been applied to model AUG ramp-down scenarios, including the impact of plasma current, auxiliary heating and shaping.Successful validation on present-day machines increases confidence that the same models, coupled to conservative assumptions regarding heat and particle transport, can be used to study optimized ramp-down strategies for DEMO.

Challenges for DEMO ramp-down scenarios
Even though the European DEMO strategy aims to maximize its reliance on a conservative physics basis that can be explored on ITER, some fundamentally new challenges will arise.
• Seeded impurities like xenon are required to radiate sufficient heat from the plasma core, maintaining the heat load to the divertor tiles manageable (even in the presence of detachment).Depending on the plasma temperature, the plasma can (locally) be in a regime where a decrease in T e leads to an increase in radiated power (from both Xe and W [14,15]), triggering a radiative collapse of the plasma.In the DEMO power balance, both the main source term (fusion power) and the main sink term (impurity radiation) are nonlinearly dependent on the plasma temperature and density, making the plasma largely self-regulated.Thus, the dynamics of a burning plasma with high radiation fraction is highly non-linear, while external actuators are less effective with respect to present-day machines.• To maintain the integrity of the thin metal wall (which must be thinner compared to ITER to allow for tritium breeding [16]), a loss of plasma control at high plasma current I p >∼ 5MA and high stored energy is unacceptable.Developing a reliable ramp-down strategy, both for planned and emergency termination of the plasma, is hence critical for the DEMO mission.
An emergency shutdown scenario in the event of divertor reattachment is discussed in [14].Divertor sweeping is proposed to delay the heat flux to the coolant becoming critical.While maintaining this emergency measure, temporarily averting target plate damage [17], a fast plasma current ramp-down is paramount.Furthermore, a fast reduction of plasma current is beneficial to reduce the forces experienced by the vacuum vessel in case of a disruption, which are proportional to I 2 p .However, note that for non-emergency ramp-up and rampdown scenarios, slower ramp phases with a slow density evolution might be preferred to allow the turbine to follow slow changes in fusion power, to maximize exploitation for electric energy production [14,18].
In this paper, the feasibility of different plasma current ramp-down rates is investigated.Through action of the central solenoid, the loop voltage at the edge of the plasma is controlled to maintain the imposed I p time evolution.The reduced loop voltage at the edge then provides a driving force for outward current diffusion.However, due to the high temperatures and the large size of a DEMO plasma, current diffusion is extremely slow.A fast current ramp-down will hence tend to peak the current density profile, or equivalently, increase the plasma internal inductance ℓ i3 , resulting in reduced controllability of the vertical position of the plasma.Since a loss of position control of the plasma column needs to be avoided throughout the entire ramp-down, the minimum time window required to safely terminate the discharge is constrained by the vertical stability limit.For the work presented in this paper, CREATE-NL simulations [19] of the vertical position stabilization control loop for DEMO [20] are used to obtain the upper limit for the internal inductance.
Note that vertical stabilization of the plasma column is projected to become more challenging for future tokamak reactors with respect to present-day devices.Future tokamaks like DEMO aim to maximize performance by significantly elongating the plasma.Compared to present devices, the current diffusion time scale is very long, while conductive walls are further away from the plasma due to the presence of tritium breeding blankets.Furthermore, measurements are complicated by the presence of 14.1 MeV neutrons and control will be less efficient due to the difficulty of putting internal coils inside the vessel and the comparatively long distance to the poloidal field coils (which are located outside of the toroidal field coils [21]).
In addition to the vertical stability limit, a set of further constraints limits the operation space of feasible ramp-downs.
• In order to ensure stable radial position control, the maximum rates of change of poloidal field coil currents impose an upper limit on the time derivative of the vertical magnetic field B v , which can be written as: Rapid changes of any of the parameters in equation ( 1), e.g. during the HL transition, can hence potentially cause a loss of radial position control.• Since the Greenwald density limit [22], the upper limit on the plasma (edge) density [23,24], is proportional to the plasma current, the density has to be decreased throughout the rampdown phase.In the present work, we assume that a constant Greenwald fraction can be maintained during the H-mode and L-mode phases.In practice, the particle confinement time and pumping capacity limit the maximum achievable density decay rate.On present-day devices, maintaining regular ELMs during the ramp-down H-mode phase is observed to be critical in order to avoid an increase in Greenwald fraction [11,25].• While terminating a burning plasma, the fusion power is reduced, by changing the isotope DT concentration and by bringing down the density (P fus ∼ n 2 e ).The presence of the inherent W impurity and the seeded Xe impurity to boost core radiation make the plasma prone to a radiative collapse: while the alpha heating drops, the average cooling factor of W and Xe increases for decreasing T e .Methods for the removal of Xe and W or large auxiliary heating power resources are required to maintain a positive power balance throughout the entire ramp-down phase.
Both the vertical and the radial position control problem illustrate the impact of the time evolution of internal plasma profiles on the magnetic control, through parameters like ℓ i3 and β pol .Conversely, the plasma shape evolution can be used as an actuator to drive changes to the plasma profile evolution: in [1], it was found that a fast decrease in plasma elongation allows to limit the increase of the internal inductance (while simultaneously widening the margin for vertical controllability).Furthermore, the plasma shape impacts the thermal confinement quality of the plasma.These examples illustrate the inherently coupled nature of the kinetic (q, T e , n e ) and magnetic (position and shape) control problems.The constraints mentioned in this section, non-linearly dependent on the plasma state itself, have to be simultaneously met.A fast transport solver like RAPTOR captures some of these nonlinearities and can hence assist in the design process of safe ramp-down strategies, as will be presented in the remainder of this paper.
The set-up of the RAPTOR simulations is described in section 3, highlighting the various non-linearities captured by the model.In section 4, a stationary operating point for a DEMO reactor is established, which will serve as the initial state for the ramp-down simulations.The time traces of auxiliary heating and Xe impurity concentration are manually optimized to find a feasible ramp-down scenario, avoiding a radiative collapse, in section 5.The upper limit for the internal inductance, obtained from CREATE-NL vertical stability calculations, is introduced in section 6.For various ramp-down rates, L-mode confinement quality assumptions and HL transition timings, the feasibility with respect to vertical and radial position control is assessed.Optimization with respect to vertical stability is studied in section 7. The RAPTOR non-linear optimization algorithms are used to optimize the ramp-down time traces of plasma current and elongation to ensure operation within the vertically stable operating envelope extracted from CREATE-NL calculations, while avoiding decreasing values of q 95 that could compromise MHD stability.The main insights resulting from our ramp-down simulations are summarized in section 8. Conclusions are formulated in section 9.

Simulation set-up for DEMO simulations
The termination simulations in this paper cover a plasma current ramp from I p flat top = 17.75MA down to I p final = 5.00MA.Depending on the I p ramp-down rate, different simulation time windows result, with t final = I p final −I p flat top dIp/dt (t 0 = 0 s at the end of flat-top).

Stationary state (initial state ramp-down)
To obtain a stationary solution, describing the plasma profiles by the end of the flat-top phase, the RAPTOR stationary state solver, described in [26], is used.The obtained stationary state, with a radially flat loop voltage profile U pl (ρ), is used as the initial state x 0 for the ramp-down simulations.

Transport equations and heat sources
• The transport equations solved for are electron heat transport T e (ρ, t), electron density transport n e (ρ, t) and poloidal flux diffusion ψ(ρ, t) (equivalent to current diffusion).These equations are evolved from t 0 = 0 s to t final on the full radial domain ρ ∈ [0 1].The boundary conditions for T e and n e at ρ = 1 are set to respectively 200 eV and 0.5 × 10 19 m −3 . 8The temperature for the ion species is set equal to the electron temperature, which is justified by the assumption that equipartitioning occurs on a fast time scale as compared to the large confinement time expected for DEMO.• For the main ions, a 50/50 deuterium/tritium fuel mix is assumed.A set of three impurities is assumed, as discussed later, with an impurity density set proportional to n e , with a user-defined, time-dependent factor.n D + n T and Z eff are solved for by imposing quasi-neutrality and evaluating the effective charge equation for Z eff .• The flat-top neutral beam deposition profiles p nbi,e and j nbi are taken from a METIS [9] DEMO simulation, with a timedependent factor scaling the profiles to match the requested total power evolution P nbi (t).The time trace of this factor is calculated before the RAPTOR simulation, taking into account the programmed plasma volume evolution.Note that while the total injected power is reduced during ramp-down, the deposition profile is maintained self-similar throughout the simulation.Updating the NBI source profiles, self-consistently with the changing plasma density and equilibrium, e.g. through direct coupling with the RABBIT code [27], is left for future work.• The EC heating is deposited in the center ρ dep = 0 with w dep = 0.1, without any current drive.• The fueling of the plasma is modeled indirectly: while the electron diffusion coefficient is set in relation to the electron heat diffusivity D e = 0.2χ e , the density gradient in the edge region of the plasma evolves in response to a dedicated term in the formula for the electron pinch velocity (see the term including µ ne in equation 2 in the Part A paper [13]).The trace of the edge density gradient is feedback controlled to track the desired time evolution for the line average density, as described in more detail in section 3.3.• The alpha power density is evaluated with the formula described in [28], consistently calculating the electron and ion heating contributions, according to the formula derived in [29].Since this alpha power model is relatively simplistic (not taking into account fast ion losses due to orbit and ripple effects and diffusion across the plasma during the slowing down time), a multiplicative factor c α = 0.73 was introduced to obtain the DEMO RAPTOR simulations reported in [30], benchmarked against more complete ASTRA simulations.
• The transition from H-mode to L-mode confinement is modeled by a gradual ramp (during the HL transition interval) of parameters in the analytical transport coefficient formulas for the T e and the n e equation (equations ( 1) and ( 2) in the Part A paper [13]).As explained in the following section 3.3, dedicated parameters enable the modification of core and edge gradients, corresponding to the confinement quality degradation and the change in core logarithmic gradients characteristic to the HL transition.

Heat and density gradient-based transport model
The present paper applies the same gradient-based transport model that has been applied for AUG ramp-down modeling in the Part A paper [13].We repeat here the main assumptions underlying the tranport model, as well as those settings that are specific to the application for DEMO.For more details regarding the analytical formulas for heat and particle diffusivities, we refer the reader to section 2.1 of the Part A paper [13].
The gradient-based model assumes the formation of three radially separated regions: (i) a central region ρ < ρ inv (= ρ q=1 ) with high transport, to mimic the profile flattening caused by sawteeth or other transport-enhancing phenomena (in the absence of a q = 1 surface, we put ρ inv = 0.1, as flattened profiles toward the magnetic axis can occur even in the absence of sawtooth activity, e.g. in the presence of kinetic ballooning modes [31]); (ii) an intermediate stiff core region characterized by constant logarithmic gradients λ Te = − d log Te dρ and λ ne = − d log ne dρ ; (iii) a pedestal region with linear gradients µ Te = − dTe dρ and µ ne = − dne dρ .
For the DEMO ramp-down simulations in the present paper, we use the gradient-based transport model because its prediction is well-defined, based on the H 98,y2 confinement factor (H 98,y2 = 1 is assumed here in H-mode, allowing for a conservative prediction), and it models across the whole plasma radius for both L-and H-mode plasmas.The LH transition is modeled through the user-defined time trace of the confinement factor and the resulting modifications of the transport coefficients.The values used for the logarithmic gradients λ Te,ne are inspired by the available literature on predictive DEMO modeling, as described later in this paper.
The time traces µ Te,ne (t), governing the pedestal confinement quality, are set as the sum of a feedforward and a feedback contribution.The feedforward contribution provides an initial estimate of the time evolution of the pedestal gradient.The feedback controller, with a proportional and an integral term, adds a corrective term to bring the plasma confinement time and line-averaged density towards pre-defined reference traces τ E ∼ H ref τ E scl and n el ∼ n el ref , making use of the error terms defined in equation (2).More technical details on application of the gradient-based transport model can be found in section 5.2.4 of [32] µ X (t) = µ ff X (t) + g p e (t) In the present work, the line-averaged density reference n el ref is set proportional to the Greenwald density n e,Gw = I p /(π a 2 ), with a user-defined factor f Gw = n el ref /n e,Gw that can be timedependent.Different confinement time scaling τ E scl could be applied; in this study the IPB98(y, 2) scaling law is used [33].
Note that the gradient-based transport model omits the need to provide a boundary condition at the pedestal top.This approach differs from the RAPTOR simulations for ITER in [26], where pedestal pressure values consistent with the EPED1-SOLPS fit reported in [34] were chosen, and from the RAPTOR simulations for AUG in [35], where pedestal pressure values were calculated with a scaling law derived based on previous experiments.In [13], the gradient-based transport model, with global constraints on the confinement scaling factor and the line-averaged density, has been successfully applied to model ASDEX Upgrade ramp-downs.
It is important to note that the IPB98(y, 2) scaling law is derived based on the data available on present-day machines, with modest levels of radiated power from the core plasma compared to a DEMO plasma.This raises the question how radiated power should be correctly accounted for in the calculation of confinement time and in the evaluation of the scaling law.In this paper, the confinement factor is calculated without subtracting the power radiated in the core.In appendix, simulations are presented to quantify the change of expected electron temperature profile and stored energy when the subtraction of core radiated power is included.

Impurity concentrations and line radiation
Three impurity species are considered: (i) an intrinsic tungsten influx is assumed since W is envisioned as plasma-facing component armour material, (ii) xenon is seeded in the core to enhance the radiated power, limiting the divertor heat load, (iii) as alpha particles born from fusion reactions thermalise, heating the plasma, they constitute a source of helium impurities.
At present, RAPTOR does not solve for impurity transport 9 .Within the simulation, the radial distribution of the impurities is taken proportionally to the electron density n e , with a user-defined time-dependent concentration factor (the resulting time evolution of Z eff is calculated self-consistently).Impurity radiation of the three impurity species is evaluated with the formula equation (3), with L imp (T e ) the impurity cooling factor taken from the ADAS database [15,37,38] P rad = n e n imp L imp (T e ) . ( Note that impurities affect the plasma power balance both by dilution of the main ions, impacting the fusion power, and by the emitted line radiation.The second process is highly non-linear with respect to temperature since the average cooling factor increases throughout the plasma core for decreasing T e (as illustrated in figure 1).Under the modeling assumption T i /T e = 1, the fusion power, and the alpha self-heating of the plasma, will decrease simultaneously.This dynamical process clearly has the potential of triggering a runaway process, with increased radiation and decreased fusion power further decreasing T e .Note that a deviation of T i /T e = 1 during ramp-down would impact the remaining alpha power, directly influencing the heat balance and margin to radiative collapse.

Equilibrium geometry
The equilibrium geometry used for the RAPTOR simulations is based on the free boundary equilibrium calculations in Plasma boundary shapes at different values of the plasma current during the ramp-down phase, simulated with the free boundary equilibrium solver CREATE-NL [19], as reported in [20].
The equilibrium geometry of these equilibria is used for the RAPTOR simulations in this paper.
CREATE-NL.The plasma boundary shapes at different times in the plasma ramp-down, at different values of the plasma current, are shown in figure 2. The geometry metrics corresponding to these equilibria are assigned to the time in the RAPTOR simulation when the corresponding plasma current is reached.For intermediate times, the geometry metrics are interpolated linearly.The CREATE-NL calculations have also been used to obtain an operating envelope for vertically stable operation during the DEMO ramp-down phase, as explained in section 6.2.
The equilibria have a lower single null configuration.The elongation is reduced during the ramp-down, while the last closed flux surface (LCFS) shape close to the X-point remains mostly unchanged, easing the heat exhaust challenge for the divertor by maintaining the magnetic geometry in the strike points region.
Note that the RAPTOR diffusion equations, as reported in [39], allow to consistently include the impact of the time derivative Φb (toroidal flux enclosed by the LCFS).Results in the present paper are slightly different from what we reported in [32], since the Φb term has been included in the present work and can have significant effects on current density peaking when the shape is modified relatively fast.

Stationary DEMO operating point
The stationary operating point established here is mainly based on considerations reported in [8,14,40].The operating parameters are summarized in table 1.A central feature of a burning DEMO plasma is the high degree of self-regulation of the plasma profiles: the power balance is dominated by the plasma self-heating by the fusion-born alpha fast particles, dependent on T i (ρ), n D,T (ρ) and fuel dilution, and the radiated power from heavy impurities, both intrinsic (W) and seeded (Xe), with a non-linear dependence on T e .Reliance on auxiliary current drive to tailor the q profile is minimized to maintain a high fusion gain Q.
• In [8], the physics-based transport model TGLF is used to predict critical temperature and density gradients.Based on those values, we set λ Te = 2 (R/L Te ∼ 6) and λ ne = 0.67 (R/L ne ∼ 2).In the L-mode phase, discussed in the next section, we assume λ Te L mode = (3/2.3)λTe H mode and λ ne L mode = (1/0.5)λne H mode , applying the same factors λ Te,ne L mode /λ Te,ne H mode as were obtained for JET and AUG in [1].The ion temperature is set T i = T e .Even though electron heating is dominant for the simulated DEMO plasma, T e ∼ T i is assumed due to the high confinement time scale with respect to the equipartition time scale in a DEMO device [41,42].• The W concentration is set to 3 10 −5 , like in [8].The Xe concentration is set to 5 10 −4 , allowing for a total radiated power of 242 MW, which is about 54% of the total heating power.This allows to limit P sep /R 0 to 22.9 MW m −1 (close to P sep /R 0 = 18.9 MW m −1 in the EU-DEMO 2018 design point reported in [14]), while a margin P sep − P LH ∼ 80 MW is maintained, to achieve good confinement and avoid any unwanted HL back transition.The helium concentration is set to 0.05, inspired by the COREDIV [43] simulation results in [8], and below the maximum He concentration of 0.075 mentioned in [30].• The fusion power is P fus = 1732 MW, below the 2 GW of the EU-DEMO 2018 design point.Note that this value is sensitive to assumptions regarding temperature and density peaking (higher reactivity in the center), Greenwald fraction, H 98y,2 factor and impurity concentrations (diluting the DT fuel).The impurity concentrations will also impact the discharge duration through the impact of Z eff on the neoclassical conductivity, and hence on the loop voltage to be provided by the central solenoid to sustain the required ohmic plasma current.• For the applied R/L ne , a density peaking n e0 /⟨n e ⟩ vol = 1.45 results, close to what has been reported in [42].With this density peaking a Greenwald fraction ⟨n e ⟩ line /n e Gw = 1.2 can be achieved, while maintaining the pedestal density about 5% below the Greenwald density (the density constraint is assumed to be active at the pedestal top location [40]).• It is interesting to note that RAPTOR predicts that bootstrap current is about half of the total plasma current.As illustrated in figure 3, the bootstrap current density is a significant source of off-axis current density (with a notable peak in the pedestal region).As a consequence, q min is relatively close to unity, with a region of low magnetic shear extending to ρ ∼ 0.25.

Ramp-down simulation with optimized heating power: discussion of actuators and constraints
In this section, we explore manually optimized time traces of auxiliary heating and the Xe impurity concentration, aiming for a feasible ramp-down scenario, avoiding a radiative collapse.Optimization with respect to vertical stability is presently not considered, but will be addressed in section 7.
The plasma current is reduced with a constant ramp rate dI p /dt = −100 kA s −1 .A time trace of the evolution of H 98,y2 is pre-defined, setting both the timing of the HL transition and assumptions regarding the confinement quality during Hand L-mode.For the simulation presented in this section, the HL transition is initialized at t = 0.2t final .The H factor transitions linearly from H 98y,2 = 1 to H 98y,2 = 0.5, over a duration 10 ∆t duration HL = 15 s.Note that these choices are relatively arbitrary at this stage.In the next section we will however perform a sensitivity study, changing the HL timing and L-mode confinement factor.An overview of the evolution of various parameters during the ramp-down simulation is given in figure 4.

Power balance and impurity concentrations
We list some of the constraints and considerations that have been taken into account when designing these ramp-down traces.

Significant auxiliary heating required throughout the entire ramp-down phase.
Even if the Xe concentration can be efficiently reduced during the ramp-down, significant auxiliary heating of the plasma throughout the entire modeled ramp-down phase (i.e. the diverted phase down to I p = 5 MA) is mandatory to avoid a radiative collapse.This is due to the combined effect of an increasing cooling factor for W and Xe for reducing T e , and a simultaneous sharp decrease of the alpha power for reducing T i and n e (the density has to be reduced simultaneously with the plasma current to avoid an increasing violation of the Greenwald density limit).Note that higher fidelity simulations are required to assess the feasibility of the proposed ramp-down rates of density and Xe impurity concentration, taking into account the expected particle confinement time (usually the particle confinement time is 5-10 times larger compared to the energy confinement time [44]), impurity confinement time and pump efficiency.
To conclude: the plasma needs significant auxiliary heating while being terminated (P ec = 50 MW+P nb = 50 MW at the beginning of ramp-down, maintaining P ec = 20 MW+P nb = 20 MW by t = t final , when I p = 5 MA).The margin with respect to radiative instability can be evaluated from the bottom right plot in figure 4: the heating power to the electrons P heat e should be maintained above the radiated power P rad , throughout the ramp-down.The margin is relatively small at the end of the HL transition: with P aux ∼ 70 MW, a margin P heat e − P rad ∼ 20 MW can be maintained.Note that the heating of discharges with high radiation during ramp-down is already common practice on present-day devices, as discussed for AUG in the Part A paper [13].Evolution of the shine-through power during ramp-down and compatibility with first wall head load limits should be assessed with dedicated simulations.

Self-consistent triggering of HL transition.
While enough heating power needs to be maintained throughout the L-mode phase of the ramp-down, the reduction of heating power during the HL transition power should obviously be significant enough to actually trigger the transition to an L-mode plasma.Note that the rapid reduction of stored thermal energy provides an effective heating term − Ẇth = −dW th /dt > 0, as the stored thermal energy crosses the separatrix, contributing to P sep = P heat − Ẇth − P rad , delaying the timing when P sep < P LH , with P LH the LH threshold power predicted by the Martin scaling law [45].
In the present simulation, self-consistency of the HL transition is ensured by making sure P sep drops below P LH during the HL transition time window between t = t HL and t = t HL + ∆t duration HL .The H factor reference trace is linearly reduced during this time window, leading to a decreasing alpha heating.Furthermore, the Greenwald fraction is changed linearly from f Gw = 1.2 to 1 during the HL transition interval.The power conducted over the separatrix P sep comes down to P LH only by the end of the HL transition time window, due to the effective Ẇth heating term.In the simulation, reducing the Greenwald fraction during the HL transition is found to increase the margin with respect to a radiative collapse, as calculated selfconsistently from the various species densities (assuming the pre-defined impurity concentration traces).

The effect of Z eff .
Note that by ramping down the impurity concentrations of He and Xe, to avoid radiative collapse, Z eff decreases throughout the ramp-down (the He concentration is reduced, which is consistent with the reducing number of fusion reactions, reducing the source term for He ions; the W concentration is maintained constant 11 ).The resulting Z eff evolution is evaluated self-consistently in RAPTOR.A reduced Z eff leads to a reduced resistivity of the plasma, slowing down current diffusion.This raises an interesting trade-off between the margins with respect to radiative and vertical instabilities: by reducing Z eff , margin with respect to a radiative instability can be improved, at the expense of a slower current diffusion, leading to a more significant increase of ℓ i3 , making the plasma more vertically unstable.Tailoring the plasma current density to limit the increase of ℓ i3 , by optimization of the time traces of plasma current and elongation, is discussed in section 7.

Heat exhaust constraints.
Finally, the reduction of impurity concentrations should not lead to a large increase in the power crossing the separatrix, that needs to be handled by the divertor.Even though this heat load would be transient, a reduced P sep with respect to flat-top values might be required due to the difficulty of maintaining the plasma detached during the ramp-down phase, as reducing density and impurity content complicate efficient dissipation in the SOL.In the 11 For the ramp-down of the ITER Q = 10 baseline scenario, predictive simulations have been performed to assess whether the influx of tungsten into the core can be avoided [7,12].Time traces of pellet fuelling and auxiliary heating are expected to have a significant impact on W accumulation.A gradual reduction of pellet fuelling was found to be beneficial to avoid the formation of significant core density gradients and the core influx of W, provided that ELM control is maintained during ramp-down.The present paper does not include predictive simulations for the level and the radial distribution of tungsten.
present simulation P sep is maintained below 240 MW throughout the H-mode phase (compared to P sep = 205 MW during flat-top burning plasma phase).This transient increase in P sep is probably beyond what can be tolerated during ramp-down.Further modeling studies into the simultaneous reduction of alpha power and radiated power, including quantitative heat exhaust constraints, are left for future work.In figure 5, the current density j par , integrated current density I encl and the loop voltage profile U pl are shown at various times during the rampdown simulation, first during H-mode and subsequently at two times during the L-mode phase.Close to the initial state of the simulation, the current density profile is relatively broad, with a significant contribution from the bootstrap current driven within the pedestal region.The simulation starts from a stationary state with a fully relaxed current density profile, characterized by a radially flat loop voltage profile.As a consequence, the ohmic current density is self-similar with the neoclassical conductivity profile.
To follow the imposed dI p /dt, the edge loop voltage is continuously reduced (except for an increase during the HL transition), as illustrated in figure 6. 12 As the edge loop voltage is reduced, to significantly negative values, the overall edge poloidal flux difference is negative.The internal current distribution of the plasma changes with a diffusion time scale that increases with the size of the device and the plasma electron temperature, leading to large values in DEMO with respect to present-day devices.In figure 6, an estimate of the characteristic resistive time is included, by evaluating τ R = µ 0 (a/2) 2 ⟨σ neo ⟩, where a is the minor radius and ⟨σ neo ⟩ is the volume average neoclassical conductivity (including a factor accounting for the impact of trapped particles [46, 47]).Clearly, the time window of the ramp-down (∼ 100s), is small with respect to the resistive time that characterizes the time scale for current diffusion, especially during the H-mode phase (τ R ∼ 1500s).
As the ramp-down is fast with respect to the current diffusion time scale, the current density evolves into an increasingly non-equilibrated state (U pl0 -U pl,edge increases, as evident in figure 6).As the loop voltage becomes peaked, the ohmic current density (which is the dominant plasma current contribution, especially during L-mode), can have a substantially different shape with respect to the neoclassical conductivity profile (and hence T 3/2 e ).The plot showing the integrated plasma current profiles (I encl in figure 5), shows that the enclosed plasma current at small radii decreases very slowly, due to the slow outward current diffusion.As the plasma volume reduces throughout the ramp-down phase, the central current density rises monotonically throughout the ramp-down (see the time trace  for j par at ρ = 0 in figure 6).Nevertheless, the total plasma current evolution (imposed as a Neumann boundary condition for the poloidal flux diffusion equation), needs to be satisfied, leading to a negative current density in the outer plasma (see the time trace for j par at ρ = 0.8 in figure 6).The MHD stability of suchlike current density profiles has not been investigated.
As the current density evolves on a slow time scale with respect to the plasma current, the plasma current ramp-down rate dI p /dt is an effective control parameter to tailor the evolution of the internal inductance ℓ i3 , as has been discussed in [11].In section 7, this feature will be leveraged to maintain the internal inductance below an upper limit given by vertical stability calculations in CREATE-NL.
In [48], a lumped parameter model for the time evolution of the tokamak plasma current and internal inductance has been proposed, with plasma resistance, non-inductive current and boundary voltage or poloidal field coil currents as inputs.The circumstances for a correlation between dℓ i3 /dt and dI p /dt are analytically derived, supporting the use of I p as a virtual actuator to control the internal inductance.

The impact of sawteeth on current density
peaking.The peaking of the current density, as described in the previous section, leads to a continuous decrease of q 0 , even though q 95 is increasing throughout the ramp-down, as shown in figure 4. As a consequence, sawteeth instabilities will likely be triggered at some time during ramp-down.In the present section, we assess the impact of sawteeth on the time evolution of ℓ i3 .In figure 7, the early HL, cold L-mode −100 kA/s ramp-down simulation is compared to a simulation where the RAPTOR sawtooth model, first presented in [49], is applied, using the sawtooth models described in [50,51].A sawtooth crash is triggered when the magnetic shear at q = 1 exceeds a user-defined critical value s q=1,crit .In this simulation we set s q=1,crit = 0.4, inspired by the values mentioned for a burning plasma in [52].The red dash-dotted traces in figure 7 show that the first sawtooth crash is triggered only in the second half of the ramp-down phase, after the plasma transitions to L-mode.The long sawtooth period during the H-mode phase of the plasma is consistent with the fast particle stabilization of sawteeth expected in a burning plasma [50, 53] (resulting in a higher expected value of s q=1,crit compared to present-day devices).The sawtooth module explicitly models the expulsion of energy and particles from inside the q = 1 radius after a sawtooth event is triggered.In the the gradient-based transport model settings (described in section 3.3), the central region with high transport has been limited to ρ inv = 0.2 for these simulations, assuming that even during the inter-sawtooth interval, core profile flattening can be expected close to the center of the plasma (e.g.like discussed in [31]).
While sawtooth crashes prevent the continuous increase of the central current density, the broadening of the current density profile, through magnetic reconnection occurring on the Alfvén time scale, is localized in the core of the plasma.The impact on the time evolution of the internal inductance is negligible, as illustrated in figure 7. Since the effect of sawtooth is modest under the present modeling assumptions, the sawtooth model is not used in the feasibility and optimization studies discussed further in this paper.

Modeling results
With the set-up discussed in the previous section, a range of different ramp-down simulations has been performed, with each of the following I p ramp-down rates: dI p /dt = −50 kA s −1 , −100 kA s −1 , −150 kA s −1 , −200 kA s −1 .Note that dI p /dt = −200 kA was used as the starting point of this study, as it is the DEMO reference design value that was also used for the CREATE-NL free boundary equilibrium control simulation in [20].For each of these ramp-down rates, two simulations are run, differing from one another in terms of HL transition timing and the assumed L-mode confinement quality: • Early HL, cold L: the HL transition is initialized at t = 0.2t final .The H factor transitions linearly from H 98y,2 = 1 to H 98y,2 = 0.5, over a duration ∆t duration HL = 15 s.These are the same assumptions that have been applied for the simulation in section 5. • Late HL, hot L: the HL transition is initialized at t = 0.4t final .
By executing the ramp-down simulations for both of these assumptions, a case with a more significant confinement transition earlier in the discharge can be compared with a more gradual confinement transition later in the discharge (giving a rough estimate for the sensitivity to the assumed H factor trace).
The resulting RAPTOR simulations are presented in figure 8.As expected, increasing the absolute value of the I p ramp-down rate leads to an faster growth of the internal inductance ℓ i3 .Furthermore, for each of the ramp-down rates, a delayed HL transition combined with an improved L-mode confinement quality leads to a more significant growth of ℓ i3 .Improved confinement during L-mode slows down the diffusion of the central plasma current, causing more significant peaking of the current density.
Let us look in some more detail at the ℓ i3 dynamics for the slow ramp-down with dI p /dt = −50 kA s −1 , shown in blue in figure 8. Interestingly, we can observe that initially, remaining longer in H-mode leads to slightly reduced values of ℓ i3 (blue dotted line) compared to the case with early HL transition (blue solid line).Note however that the increase in growth rate of ℓ i3 after the early HL transition is very modest when compared to e.g. the dynamics observed for AUG in Part A [13].While the late HL transition, hot L-mode case initially maintains lower ℓ i3 values, the higher L-mode confinement during the second half of the ramp-down phase eventually leads to a more significant increase of ℓ i3 with respect to the early HL transition, cold L-mode case.
Furthermore, increasing the I p ramp-down rate leads to increasingly negative values of the parallel current density near the edge of the plasma, as shown in figure 9.The −200 kA s −1 ramp-down with late HL transition features the largest ratio between average resistive diffusion time τ R and the rampdown time t final , and the most negative values of j par .Looking at the profile of the enclosed plasma current at the final time step of the simulation (when I p = 5 MA), we observe a significant negative current with a magnitude of about 3 MA in the outer plasma region ρ > 0.7.As negative edge loop voltages are applied for all modeled ramp-downs (except for a small time interval in the −50 kA s −1 ramp-down with early HL transition), the expected edge poloidal flux difference is negative for all cases.Interestingly, these simulations do not indicate a significant difference of the edge poloidal flux swing that is required for the different I p ramp-down rates, as quantified in panel (d) of figure 8.Note that the edge poloidal flux difference is not equal to the flux swing of the central solenoid, as the flux due to the plasma external inductance and the flux due to changing currents in the poloidal field coils need to be accounted for (especially when changing the plasma shape).However, these simulations indicate that significant recharging of the central solenoid can be expected during ramp-down for DEMO.
To ensure radial position control, the time derivative of the vertical magnetic field must stay below an upper limit, depending on coil voltage limits imposed by power supplies and/or the superconductor.While we have presently no upper limit value available, we evaluate dB v /dt with equation (1) for the different simulations, as shown in figure 8.The time derivative dB v /dt reaches the largest absolute values during the HL transition (as expected from the β p dependence in equation ( 1)).Interestingly, the best-case scenario H 98y,2 trace regarding vertical stability (early HL, cold L) is the most demanding regarding radial position control, with the most significant peak in Table 2. From free boundary equilibrium control calculations with CREATE-NL [19], the above combinations of (Ip, κ, ℓ i3 ) are considered controllable, as reported in [20].For a given plasma current, the elongation and internal inductance values provide an upper constraint on the stable operating envelope.the dB v /dt trace.To disentangle the effects of HL timing and L-mode confinement on the dB v /dt trace, we also executed a simulation with late HL transition to cold L-mode (for the −100 kA s −1 case).For this additional simulation, dB v /dt values similar to the early HL transition to cold L-mode case have been obtained, indicating a modest impact on radial position control of the HL timing itself (under the assumption that the Greenwald fraction during H-mode can be maintained constant).

Comparison of RAPTOR-predicted ℓ i3
versus CREATE-NL vertical stability limit 6.2.1.CREATE-NL vertical stability limit.In [20], free boundary equilibrium control calculations with CREATE-NL [19] are presented for the diverted phase of the DEMO plasma ramp-up and ramp-down.For the ramp-down, the limits to a vertically stable operating envelope are mapped out with respect to the plasma current I p , elongation κ and internal inductance ℓ i3 .The resulting sequence of triplets (I p , κ, ℓ i3 ) of controllable operating parameters is repeated in table 2. Each of these three parameters play an important role in the assessment of vertical stability: while a degraded vertical control efficiency can be anticipated for an increased internal inductance ℓ i3 , corresponding to a more significant peaking of the current density, this tendency can be counteracted by adjusting the plasma shape, reducing the elongation.Furthermore, controllability of the vertical position improves at lower plasma current, as the coil currents for vertical stability control in the CREATE-NL model for DEMO are more effective to counteract vertical position excursions at lower plasma current [20].
Considering the upper limits for internal inductance ℓ i3 and elongation κ for a given plasma current I p (as summarized in table 2), the potential of different ramp-down rates to maintain the internal inductance below the upper limit for vertical controllability can now be assessed.For the purpose of the feasibility study presented here and the optimization introduced in the next section, we employ these values to extract a constraint on the internal inductance, dependent on the plasma current value.We assume that the plasma can be maintained vertically stable if ℓ i3 < f margin ℓ i3 CREATE (I p ).We include a margin factor f margin > 1, as the present constraint ℓ i3 CREATE (I p ) is a conservative assessment, not an optimized limit.Various considerations could increase the maximum allowed ℓ i3 .
• Optimizing position and shape control, while using the consistent kinetic profile evolution, is expected to provide some improvement with respect to vertical controllability.• Inclusion of in-vessel coils in the DEMO design would facilitate more effective vertical control.• A faster decrease of elongation κ with respect to the reference in [20] improves vertical controllability for a given plasma current.
The vertical stability constraint ℓ i3 CREATE (I p ) is indicated in figure 10, with the diamond symbols representing the equilibria from CREATE-NL that are introduced in table 2. Allowing for some margin with respect to this conservative CREATE-NL result, the lines with margin factor f margin = 1.15 and 1.30 are also shown.The dark grey region with ℓ i3 < ℓ i3 CREATE (I p ) is the region where vertical stability can be guaranteed, while the white region ℓ i3 > 1.30ℓ i3 CREATE (I p ) contains operating points that can likely not be maintained vertically stable.Note that the x-axis is −I p , so that moving to the right on the abscissa corresponds to progressing time during the ramp-down phase.Later in the ramp-down, at lower plasma currents, a larger internal inductance can be maintained vertically stable.

Comparison to modeled ℓ i3 traces.
Let us now superimpose the ℓ i3 traces modeled in RAPTOR on the  3.As a first conclusion, we can observe that the tendency of the plasma to peak the current density during the ramp-down phase leads to a more significant increase of ℓ i3 compared to the increase of the constraining value ℓ i3 CREATE , even for a conservative ramprate of dI p /dt = −50 kA s −1 .However, allowing for some margin on the vertical stability constraint (justified by the reasons listed earlier), the dI p /dt = −50 kA s −1 ramp-down simulations stay below upper limits with f margin respectively 1.15 and 1.30.For the fastest ramp-down assumption (dI p /dt = −200 kA s −1 ), both H factor trace assumptions lead to a severe violation of the upper ℓ i3 limit with f margin = 1.30.Also for the simulations with intermediate I p ramp-down rates (dI p /dt = −150 kA s −1 , −100 kA s −1 ), the upper ℓ i3 limit with f margin = 1.30 is violated.In the following section, we will attempt to optimize the plasma current trace and the plasma shaping evolution to bring the ℓ i3 trace below the upper constraint, with f margin either 1.15 and 1.30.

Optimized I p and shaping evolution to avoid vertical instability
In the present section, feasible DEMO ramp-down scenarios are searched for by numerical solution of an optimal control problem for the ramp-down phase.The mathematical formulation and solution procedure is discussed in section 7.1, while the actual results are presented in section 7.2.Loosely speaking, the aim is to ramp down the plasma current I p as quickly as possible, while satisfying all physical and technical constraints.Since the plasma density should decrease proportional to I p to avoid density limit excursions and since the fusion power is proportional to the square of the plasma density, a fast decrease of I p corresponds to a fast decrease of the plasma thermal energy (also the confinement time is expected to decrease proportional to I p ).A fast I p decrease hence reduces the potential impact of depositing the stored thermal energy on the reactor vessel first wall.Furthermore, the electromagnetic forces acting on the vessel after plasma disruption are proportional to I 2 p .To reflect these considerations, the time integral of the plasma current is chosen as a cost function: The objective of the optimization problem is to minimize J Ip .
The factor ν Ip is chosen to have a cost function around unity for the initial condition for the ramp-down traces.

Constraint functions.
An extensive set of physical and technical limits constrain the ramp-down phase: • Density limit: as the time evolution of the plasma current I p (t) is adjusted during the optimization, the reference density trace n el, ref (t) (defined in section 3.3) is updated after each iteration step, to keep the allowed Greenwald fraction time trace unchanged throughout the optimization procedure.The present work does not assess whether sufficient pumping capacity is available to reduce the density at the required rates.• Vertical stability limit: ℓ i3 < f margin ℓ i3, CREATE (I p ), as discussed before.Note that by optimizing the plasma current, both the internal inductance ℓ i3 and the upper limit ℓ i3, CREATE (I p ) can be adjusted.The corresponding constraint in the optimization problem is formulated as c 1 = ℓ i3 − f margin ℓ i3 CREATE (I p ) < 0 and should be enforced over the entire ramp down interval.• Ideal MHD stability: to maintain a margin with respect to ideal MHD limits, a constraint is added to enforce a monotonic increase of q 95 while q 95 < 4.5 (assuming that when q 95 > 4.5, a decrease can be safely allowed).This constraint is included as plasmas with a lower q 95 value are more prone to MHD instabilities, since they operate closer to the ideal MHD limit.Note that this is not a hard limit and could be relaxed.The corresponding constraint in the optimization problem is formulated as c 2 = − q95 < 0 and should be active only for those time points when q 95 < 4.5.
Through the formalism introduced in [54], the state constraints c 1 (x(t)) < 0 and c 2 (x(t)) < 0 (dependent on the plasma state x(t), which contains the radial distribution of T e , n e and ψ, as evolved by RAPTOR) are formulated as integral constraints (instead of a separate constraint for each t k ): The weight w i is set to unity for those times where the corresponding constraint is active (for those times when q 95 < 4.5 for c 2 ).We obtain the two integral constraints: C ℓi3<f margin ℓi3 CREATE(Ip) ⩽ 0 and C q95 >0 ⩽ 0.

Optimization variables.
The actuator traces that are optimized to minimize the cost function, while satisfying the constraints, are the plasma current I p (t) and the elongation of the plasma κ(t).Within the optimization problem formulation, the actuator traces u i (t) (u 1 (t) = I p (t) and u 2 (t) = κ(t)) are parametrized by the vector [p 1,1 . . .p 1,n1 p 2,1 . . .p 2,n2 ] T .By multiplication with a set of piecewise linear basis functions (respectively n 1 functions P 1j (t) and n 2 functions P 2j (t)), followed by summation, the actuator traces u i (t) are recovered: The initial and final values of I p and κ are maintained unchanged by fixing the parameters p 1,1 , p The timing of the HL transition is maintained fixed throughout the optimization.The simulations in section 6 and AUG experiments in Part A [13], as well as cross-machine analysis in [11], highlight the importance of this parameter in determining the ramp-down dynamics.In the present work, sensitivity to this parameter has been tested by running optimizations for different assumptions regarding the HL timing.A natural follow-up of the present work would be the inclusion of the HL timing as an optimization variable, as demonstrated in section 4.4 of [55].Auxiliary current drive, e.g.ECCD, has not been considered as an actuator to optimize the current density profile in the present work.Note that EC deposition control could also be needed to avoid tungsten accumulation in the core, or for NTM control [10].

Using the elongation trace κ as an optimization variable.
As the equilibrium impacts the radial transport equations through geometric coefficients in the diffusion PDEs, the time evolution of the elongation κ can be used as an optimization variable, to impact the plasma state dynamics.In RAPTOR, these metrics, stacked together in the vector g, are calculated from the output files of a CHEASE equilibrium solution.From the CREATE-NL equilibrium solutions introduced in section 3.5, a sequence of metrics g κ for seven different values of plasma elongation κ is obtained.During a ramp-down with constant ramp rate dI p /dt, g for intermediate time points is evaluated through linear interpolation in time.
We attempt to optimize the time evolution of the plasma shape, by finding the optimum time trace of the elongation κ(t).The fact that the dependencies of the geometric factors g on κ are not analytically available, has consequences for the optimization routine: • Cost and constraint function gradients ∂J ∂p and ∂C ∂p have to be evaluated numerically with finite differencing and cannot be evaluated analytically, slowing down the optimization.
• Lacking a coupled Grad-Shafranov solver, a consistent evaluation of the updated metrics g for a changed κ cannot be obtained within the optimization routine.
After each iteration step, the geometric factors g have to be recalculated to match the new trial of the time evolution of the elongation κ.A linear interpolation scheme is applied to update the metrics g.The two reference equilibria with most similar elongation κ are identified (out of the set of CREATE-NL equilibria κ ref ,i ): Then, a linear interpolation in κ is applied to obtain the metric vector g: Once an optimized ramp-down scenario is found, the adequacy of this linear interpolation scheme can be checked by running a set of CHEASE simulations for time points along the optimized trajectory, as illustrated in section 7.4.
7.1.5.Ramp-down optimal control problem.To summarize, the ramp-down optimal control problem can be written as: A ineq p ⩽ b ineq (actuator limits) (8d) C ℓi3<f margin ℓi3 CREATE(Ip) ⩽ 0 and C q95 >0 ⩽ 0 (state constraints).(8e) The RAPTOR state evolution equation f defines the time evolution of the plasma state x(t).The actuator limits equation (8d) are included to impose both I p (t) and κ(t) to be monotonically decreasing.Note that this condition can be written as a linear inequality constraint on p, as du i /dt ⩽ 0 translates directly to ∑ ni j dP ij (t)/dtp i,j ⩽ 0. At this stage, no minimum time derivative is set.Position and shape control studies including the poloidal field coil currents are required to establish how fast plasma current and plasma shaping can be changed during the ramp-down phase.
The algorithm applied to solve the non-linear, constrained optimization problem formulated in equations (8a)-(8e) is sequential quadratic programming (SQP) [56], as implemented in the Matlab function fmincon.This algorithm was applied before in [1,54].Even though cost and constraint function gradient information has to be obtained numerically, the fast run time of a single RAPTOR simulation allows to maintain the full solution of the non-linear optimization problem computationally tractable (a few hours on a single CPU).
Generally speaking there is no guarantee that the obtained optimum trajectory is a global optimum.Rerunning the optimization from different initial conditions allows to increase confidence in the obtained optimum solution.Finally, the selection of the number of optimization variables needs to give the optimizer enough degrees of freedom to be able to find an optimum, while a too high dimensionality leads to the risk of over-fitting with respect to specific settings of model parameters.Verifying the obtained optimum trajectory with more complete integrated modeling tools could increase confidence in the obtained optimum.

Optimized DEMO ramp-down scenarios
A set of optimized DEMO ramp-down scenarios is summarized in figure 11, with different assumptions regarding the total ramp-down time window and the time evolution of the reference H 98y,2 : (i) dI p /dt = −100 kA s −1 ; late HL transition, hot L-mode; (ii) dI p /dt = −150 kA s −1 ; early HL transition, cold L-mode; (iii) dI p /dt = −200 kA s −1 ; early HL transition, cold L-mode.
The respective trajectories given as initial conditions to the optimizer are the corresponding RAPTOR simulations shown in figures 8 and 10.The red dash-dotted traces in figure 11 represent the optimum found for the optimization problem formulated in equations (8a)-(8e), with C ℓi3<f margin ℓi3 CREATE(Ip) ⩽ 0 and C q95 >0 ⩽ 0 as active constraints.The second constraint is added as a significant reduction of the elongation during the early ramp-down can lead to low q 95 values that could compromise ideal MHD stability.To counteract this effect, dq 95 /dt > 0 when q 95 < 4.5 is imposed, ensuring q 95 increases monotonically during the early ramp-down.
For cases (1) and ( 2), a feasible solution is found with f margin = 1.15, while for case (3), f margin has to be raised to 1.30 to find a feasible solution.Both I p and κ are approximated by linear segments, parametrized by three intermediate knot points (at the intersection of linear segments).The end points are maintained equal to the initial condition, while the values of the knot points constitute the optimization variables (6 optimization variables in total).
For each of the three cases, the optimized I p (t) and κ(t) trajectories feature both a sharp initial decrease in the first segment of the ramp-down.The initial fast reduction of I p leads to a very fast increase of ℓ i3 .However, the fast decrease of I p also leads to a fast increase of the upper limit f margin ℓ i3 CREATE (I p ).Furthermore, the decrease of the I p and κ ramp-rates at the first knot point leads to a knee point in the ℓ i3 (t) trace.For the three cases, an optimum trajectory is found that avoids a violation of the vertical stability constraint throughout the entire rampdown phase, while featuring a monotonic increase of q 95 for q 95 < 4.5.
To illustrate the importance of the C q95 >0 ⩽ 0 constraint, the green dashed traces in figure 11 represent the ramp-down traces that are obtained if the optimum I p trace is applied, while reducing κ immediately to the minimum value during the first time segment.A fast reduction of the elongation is beneficial to increase the vertical controllability of the plasma.Additionally, we observe that the fast κ reduction leads to a significant reduction of the internal inductance, providing a larger margin with respect to the upper limit f margin ℓ i3 CREATE (I p ).However, the fast reduction of elongation leads to a significant decrease of q 95 , potentially compromising MHD stability (as observed in the ASDEX Upgrade experiment described in [13]).
We conclude that the non-linear optimization routine manages to find a feasible ramp-down trace, optimizing the I p and κ traces to satisfy both constraints C ℓi3<f margin ℓi3 CREATE(Ip) ⩽ 0 and C q95 >0 ⩽ 0. In the next section, we investigate in more detail the impact of respectively I p and κ on the time evolution of the internal inductance.

Interpretation of the impact of the optimum Ip and κ traces on ℓ i3
To understand better the dynamics underlying the obtained optimum trajectory, we study in more detail the optimized I p and κ traces for the case with dI p /dt = −100 kA s −1 , late HL transition, hot L-mode in figure 11.Various simulations are performed to understand both the individual and the joint impact of the optimized I p and κ traces, as presented in figure 12 and explained below.
• Initial ramp-down trajectory with constant dI p /dt and dκ/dt (black solid trace) • Optimized shaping evolution κ(t), while maintaining the initial I p (t) evolution with constant dI p /dt (blue dotted trace) The faster reduction of elongation κ leads to a reduction of ℓ i3 compared to the initial ℓ i3 trace, throughout the entire simulation time window.Note however that the value of q 95 reduces below its flat-top value, potentially compromising ideal MHD stability.• Optimized plasma current evolution I p (t), while maintaining the initial κ(t) evolution with constant dκ/dt (red dashed trace) When maintaining the plasma shaping evolution unchanged, the fast plasma current reduction leads to larger values of of the internal inductance ℓ i3 (t) throughout most of the simulated time window, even though the final value of ℓ i3 is reduced.The reduction of the I p ramp-rate around I p = 10 MA leads to a knee point in the internal inductance ℓ i3 trace.This illustrates that I p is an effective actuator to tailor the time evolution of ℓ i3 .For the same shaping trajectory, the lower plasma currents lead to a more significant increase of q 95 .• Optimized plasma current I p (t) and shaping evolution κ(t) (green dash-dotted trace) Application of the individual optimized traces of I p and κ is insufficient to maintain ℓ i3 < 1.15ℓ i3 CREATE (I p ).However, as shown in the right panel of figure 12, the combined impact of both actuators is successful in satisfying the constraint, while q 95 is monotonically increasing, remaining above the flat-top value.

Consistency of the optimized evolution of shaping and kinetic profiles with equilibria
Once an optimized ramp-down scenario is found, a consistent set of CHEASE equilibria can be calculated for time points along the optimized trajectory.This way, one can validate whether the impact of the optimized κ trace on the plasma state evolution has been correctly captured by the interpolation scheme to update the metrics g introduced in equation ( 7).An automated function allows the user to provide the time points for which a consistent equilibrium is desired.For these times, equilibria are calculated with consistent q 95 , elongation κ and kinetic profiles (p ′ and TT ′ consistent with the pressure p and current density j par of the optimized RAPTOR simulation are provided to CHEASE).To prepare a LCFS contour with the desired value κ opt for CHEASE, a set of plasma boundary control points (R opt,i , Z opt,i ) are defined.These control points are calculated by selecting the LCFS contour (R i , Z i ) of the equilibrium from the original sequence with elongation κ closest to the desired value κ opt and by subsequently rescaling the vertical distance to the (SND) X-point for each control point: The initial sequence of CREATE-NL equilibria and a sequence of equilibria along an optimized ramp-down trajectory are indicated by sets of diamond symbols in figure 13.The boundary shape for an equilibrium at t = 55 s (identified by the green boxes in the I p and κ time trace plots) is compared between the original equilibrium (blue solid line) and the updated, consistent equilibrium (red dash-dotted line).11).The values of plasma current and elongation of the equilibria initially obtained from CREATE-NL are indicated by the blue diamonds.The red diamond symbols indicate the equilibria calculated during post-processing of the optimized results.These equilibria are calculated with values of q 95 , κ and kinetic profiles consistent with the optimal trajectory.For an equilibrium at t = 55 s (identified by the green boxes): the initial equilibrium boundary shape obtained from CREATE-NL is compared to the boundary shape with optimized κ.Now that a set of consistent equilibria has been calculated, we can verify whether the simplified treatment of equilibrium adjustments in the optimization scheme can be justified.Let us compare the ℓ i3 and q 95 traces obtained with the optimizer (applying a simple, linear interpolation technique to evaluate g) with the respective values obtained with a RAPTOR simulation using the fully consistent CHEASE results as underlying equilibria.For ℓ i3 the changes are modest, justifying the applied linear interpolation procedure.For the q 95 trace, the change between the optimized RAPTOR simulation without (black dash-dotted) and with consistent equilibria (green dotted) is more significant.The assessment of MHD stability of the proposed ramp-down trajectory requires more in-depth analysis, left for future work.

Summary of the main findings for DEMO ramp-down optimization
Let us summarize the main findings we can extract from the optimization of the diverted phase of the DEMO ramp-down in this paper, and from the insights gained in Part A [13]: Reducing the plasma current ramp-rate |dI p /dt| • (+) slows down the increase of ℓ i3 , increasing the vertical stability margin (this effect is more efficient in conjunction with reducing the plasma cross-section); • (+) reduces the magnitude of negative edge currents, which is in general safer with respect to ideal MHD stability; • (+) increases the margin with respect to the Greenwald limit, relaxing density decay requirements; • (+) unlike AUG observations, the poloidal flux swing for DEMO is little affected by changes in the length of the rampdown time window and recharging is observed for all considered cases; • (−) prolongs the ramp-down duration and the dwell time at high I p .

Delaying the HL transition
• may be necessary if the fusion power cannot be reduced sufficiently rapid; • (+) reduces flux consumption (more recharging), although the difference becomes smaller for faster ramp-down rates; • (−) demands solutions for density decay and impurity control (fueling, ELM control . ..); • (−) demands a significant amount of installed auxiliary power, to maintain a sufficient P sep as P alpha reduces; • (−) prolongs the dwell time at high W th .

Maintaining L-mode heating
• (+) helps maintaining a positive power balance, increasing the radiative collapse margin after the HL transition; • (+) increases the density limit margin; • (+) reduces the drop of β pol during the HL transition, easing radial position control; • (−) while a sudden increase of ℓ i3 is averted by maintaining heating for AUG L-modes, the highly transient current diffusion dynamics in the DEMO ramp-down show a dominant impact of the loop voltage profile shape rather than the T 3/2 e profile; as increased temperatures lead to a larger resistive time τ R , slowing down the current diffusion, larger values of ℓ i3 and a larger magnitude of negative edge currents result, reducing the vertical stability margin; • (−) the reduction of separatrix power needs to be sufficiently significant to trigger the HL transition (considering the effective heating from −dW th /dt), to be managed by real-time control.

Reducing the plasma cross-section (reducing the elongation κ)
• (+) slows down the increase of ℓ i3 , increasing the vertical stability margin (this effect is more efficient in conjunction with reducing |dI p /dt|); • (+) reduces the vertical instability growth rate by reducing the elongation, increasing the vertical stability margin; • (−) reduces q 95 , reducing the MHD stability margin; • should maintain the magnetic geometry in the strike point region unchanged, as well as the shape in front of antennas if ion cyclotron (or lower hybrid) heating is being used.

Reduction of heavy impurity content of the plasma core
• (+) reduces the power radiated from the plasma core, increasing the radiative collapse margin; • (−) increases conductivity and the resistive time τ R , slowing down current diffusion, leading to larger values of ℓ i3 , reducing the vertical stability margin; • (−) increases the separatrix power to be handled by the SOL.

Conclusion
The safe termination of burning plasmas is of crucial importance for the exploitation of DEMO, as very few disruption events can be tolerated (especially at plasma currents above I p ∼ 5 MA).The high radiated power fraction of DEMO, to limit the heat load to be handled by the divertor, will make the scenario sensitive to excursions from a nominal scenario, both in stationary state and during transient phases like ramp-down: a decrease of electron temperature leads to increasing average cooling factors for the intrinsic tungsten species and the seeded xenon species, potentially triggering a runaway process towards a radiative collapse.To avoid a radiative collapse during the ramp-down phase, significant plasma heating needs to be maintained during the L-mode phase (with the level of heating depending on how efficiently impurities can be removed from the plasma).
The present paper applies the gradient-based transport model introduced in [1].A stationary RAPTOR reference for the DEMO stationary operating point is established, based on operating conditions discussed in the literature [8,14,40].The simulation includes the effect of plasma dilution by the fusion-born helium species as well the radiation from tungsten and xenon, making use of ADAS cooling factor data [37].A reactor operating point has been presented, considering the active constraints on performance, as well as physics limits.Importantly, a DEMO plasma is characterized by a high degree of self-regulation of the kinetic profiles (posing a challenge for kinetic profile control): the power balance is dominated by the plasma self-heating by the alpha fast particles, dependent on temperature and density profiles and fuel dilution, and the radiated power from heavy impurities, with a non-linear dependence on T e .
The stationary operating point is used as initial condition for a series of ramp-down simulations.The technical and physical constraints that need to be simultaneously satisfied throughout the entire ramp-down phase raise a set of tradeoffs when setting the actuator time traces.Leveraging the fast run time of the code, feasible DEMO ramp-down scenarios are developed in RAPTOR.A quantitative estimate is presented of the significant auxiliary heating required throughout the ramp-down phase, to avoid a radiative collapse.
Fast ramp-down scenarios are critical for emergency shutdown of the burning plasma, e.g. in case of divertor reattachment [14], and are actively analysed on present-day machines, e.g. after large tearing modes [2].Ramping down the plasma current tends to cause a peaking of the current density profile, with a faster plasma current ramp-down rate leading to more significant peaking.In our simulations, significant peaking of j par with inversion of the plasma current direction near the plasma edge is routinely observed.The corresponding increase of the plasma internal inductance ℓ i3 poses a challenge for the vertical position controllability of the plasma column, as thoroughly assessed in this paper.The MHD stability of these current density profiles should be assessed in future work.
An upper limit on the internal inductance from CREATE-NL [19] free boundary equilibrium control calculations is introduced [20], dependent on plasma current I p and elongation κ.The feasibility of different plasma current rampdown rates with respect to this vertical control limit is assessed: while a ramp-rate of dI p /dt = −50 kA s −1 (with ramp-down duration t final = 256 s) seems conservative, ramprates faster than dI p /dt = −100 kA s −1 (t final = 128 s) require active optimization.For a given ramp-rate, the increase of ℓ i3 can be limited by enhancing the outward diffusion of the current density: a lower temperature (lower L-mode confinement) or an increased effective charge Z eff allows for a faster current diffusion time scale.The time traces of the plasma current I p (t) and the plasma elongation κ(t) provide effective actuators to tailor the internal inductance trace ℓ i3 (t), allowing furthermore to reduce the maximum ℓ i3 without extending the time window used for the ramp-down.A reduction of the plasma elongation allows to enhance the vertical controllability of the plasma column, while a faster compression of the plasma results in a reduced increase of ℓ i3 according to RAPTOR modeling.Furthermore, during fast shape changes, a significant impact of Φb (time derivative of the toroidal flux enclosed by the LCFS) on current density peaking has been identified.However, while reducing the plasma elongation, q 95 should be maintained sufficiently high, by constraining how fast the plasma shape can be changed.
Applying the RAPTOR automated optimization routines introduced in [1], a set of feasible ramp-down trajectories has been derived, respecting dual constraints: while the internal inductance ℓ i3 is limited below the upper limit from CREATE-NL, a monotonic increase of q 95 during the early ramp-down is forced.Imposing ℓ i3 < f margin ℓ i3 CREATE (I p ) with f margin = 1.15, a feasible ramp-down strategy could be found for an average ramp-down rate of dI p /dt = −100 kA s −1 (t final = 128 s), even for conservative assumptions regarding the HL transition and confinement quality (for dI p /dt = −150 kA s −1 , t final = 85 s, a feasible ramp-down with f margin = 1.15 could only be found under favorable confinement conditions: an early HL transition and low confinement quality during L mode allowing for faster outward current diffusion).For dI p /dt = −200 kA s −1 (t final = 64 s), a feasible ramp-down scenario could be found under favourable confinement conditions and assuming a larger margin f margin = 1.30 with respect to the constraint from CREATE-NL can be achieved.Note that the avoidance of negative current densities near the plasma edge has not been included as an objective in these optimizations.While geometry modifications between iterations of the optimizer have been implemented in an ad-hoc way, this approach has been successfully validated by re-simulating the optimized trajectories with a loose coupling to the CHEASE fixed boundary equilibrium solver.
The inherently coupled nature of the kinetic (q, T e , n e ) and magnetic (position and shape) control problems demands integrated simulations to test the feasibility of established ramp-down strategies.The kinetic profile evolution of the optimum ramp-down trajectory should be confirmed with higher fidelity integrated modeling tools (including impurity transport, pedestal model etc).Ideally, whole-device modeling should be envisioned, including SOL dynamics and pumping efficiency.Iteration or coupling between a transport solver and a free boundary equilibrium control model (like CREATE-NL [19]) is required to assess controllability.More specifically: for the optimized ramp-down trajectories proposed in this paper, it should be verified whether the poloidal field coils allow for the proposed changes in plasma elongation, and whether the plasma equilibria can be maintained radially and vertically stable (for the modeled kinetic profile evolution, plasma elongation and for a given DEMO design, e.g. with or without in-vessel coils).
It should be emphasized that the optimized trajectory can be tailored in real-time with accurate knowledge, since RAPTOR is RT-compatible on present-day tokamaks and much more time and computing power is available in a DEMO reactor.In particular the best timing for the H-L transition can be determined under several pre-defined conditions, including the Greenwald fraction time evolution, and triggered in real-time when approaching an unsafe operating boundary.Similarly, the level of heating in L-mode should be as low as possible, while controlling the adequate power balance with respect to the radiated power.

Appendix. Evaluating τ E and τ E scl in a plasma with high levels of radiated power
When evaluating the confinement enhancement factor H = τ E /τ E scl , a proper definition of the loss power P L is required, present both in the confinement time formula τ E = W th /P L , and in the scaling law dependency τ E scl = τ E 98y,2 ∼ P −0.69 L .The standard implementation of the gradient-based transport model assumes P L = P heat − Ẇth , with P heat the total plasma heating power (ohmic, alpha and auxiliary) and Ẇth = dW th /dt (this term is equal to zero for the stationary operating points shown in this appendix, but non-zero during ramp-down).This definition of the loss power does not subtract any fraction of the power radiated directly from the core plasma.In [57], IPB98(y, 2) confinement scaling predictions are compared to ASTRA-TGLF simulations to come up with a proper P L correction term accounting for the core radiated power.By subtracting 60% of the radiated power inside ρ = 0.75, the best agreement between simulation and scaling law is obtained, i.e.P L = P heat − Ẇth − P rad core with P rad core = 0.6 ´0.75 ρ=0 p rad dV.Applying these corrections in the gradient-based model impacts the predicted temperatures and stored thermal energies.For a given H factor, e.g.equal to unity, one can write τ E = W th /P L ∼ τ E scl = τ E 98y,2 ∼ P −0.69 L , hence for the resulting stored energy one can write W th ∼ (P L ) τE (P L ) −0.69 τ E scl .While applying the P rad core subtraction in (P L ) τE leads to a reduced W th , the loss power correction in the scaling law tends to increase W th (since confinement degrades with increasing power).Depending on which of both loss power factors (P L ) τE and (P L ) τ E scl is corrected with a radiated power subtraction, we can hence distinguish four cases 13 (the three cases for which a RAPTOR stationary state is found are shown in figure A1 and  table A1): • The electron temperature profile for the reference case without loss power corrections is represented in blue in figure A1. • The net effect of the double loss power correction proposed in [57] leads to a reduced W th .This can be illustrated by comparing the blue and the dashed red T e profiles in figure A1.Table A1.The stationary DEMO operating points obtained in RAPTOR are presented, for different assumptions regarding the way radiated power is taken into account when calculating the loss power.The three cases have an H factor equal to unity H 98y,2 = (P L ) 0.69 τ E scl W th /(P L )τ E = 1.The first two columns indicate whether or not P rad core is subtracted, respectively in the loss power terms (P L )τ E and (P L )τ E scl . (

Figure 1 .
Figure 1.Radiated power from three impurities, for concentrations n He /ne = 0.05, n W /ne = 3 10 −5 , n Xe /ne = 5 10 −4 , as evaluated with ADAS cooling factor data, with respectively the end of flat-top Te profile (solid lines) and a L-mode Te profile (dashed lines).Xenon, the seeded core impurity, is the dominant radiator.Note how the HL transition during the ramp-down leads to a significant increase in radiated power, from both W and Xe.A combination of plasma heating and Xe removal is required to avoid a radiative collapse, as studied later in this paper.(a) Te(ρ); (b) p rad (ρ).

Figure 2 .
Figure 2.Plasma boundary shapes at different values of the plasma current during the ramp-down phase, simulated with the free boundary equilibrium solver CREATE-NL[19], as reported in[20].The equilibrium geometry of these equilibria is used for the RAPTOR simulations in this paper.

Figure 5 .
Figure 5. Illustration of current diffusion dynamics during ramp-down, highlighting radial profiles at three different times during the ramp-down.(a) Ip(t) (the three vertical lines correspond to the times for which the radial profiles are shown in (d)-(f)); (b) β N (t); (c) ℓ i3 (t); (d) parallel current density jpar(ρ) at three times during ramp-down, corresponding to Ip = 17.5MA,Ip = 10.0MA and Ip = 5.5MA; (e) idem for the enclosed current I encl (ρ); (f) idem for plasma loop voltage U pl (ρ).

Figure 7 .
Figure 7.The impact of the sawtooth model on the early HL, cold L-mode, −100 kA s −1 ramp-down simulation is illustrated.The simulation without sawteeth is shown in blue solid lines, while the simulation with sawteeth and critical magnetic shear s q=1,crit = 0.4 is shown in red dash-dotted lines.For the central temperature, a zoom panel is included to highlight the evolution of T e0 after a sawtooth crash.Only a minor impact on the ℓ i3 evolution is observed.(a) q 0 (t); (b) ℓ i3 (t); (c) radius of the q = 1 surface ρ q=1 (t); (d) Te(ρ, t) at ρ = 0 and ρ = 0.5.

Figure 10 .
Figure 10.Time traces of ℓ i3 within the (−Ip, ℓ i3 ) plane, for the RAPTOR simulations shown in figure 8.The stability limit f margin ℓ i3 CREATE (Ip) allows to assess the margin with respect to the vertical stability limit obtained in CREATE-NL.

Figure 11 .
Figure 11.Initial and optimized time traces for three DEMO ramp-downs: (1) dIp/dt = −100 kA s −1 ; late HL transition, hot L-mode; (2) dIp/dt = −100 kA s −1 ; early HL transition, cold L-mode; (3) dIp/dt = −200 kA s −1 ; early HL transition, cold L-mode.The time traces of the initial ramp-down simulation are shown in solid blue lines.The time traces of the ramp-down simulation with optimized Ip and κ trajectories are shown in red dash-dotted lines.The green dashed lines correspond to a ramp-down simulation with an Ip trace corresponding to the optimum, and a fast reduction of the elongation κ, leading to a violation of the constraint on q 95 .Feasible ramp-down traces are found with f margin = 1.15 for (1) and (2) and with f margin = 1.30 for (3).(a) Ip(t); (b) elongation κ(t); (c) q 95 (t) (the grey area indicates q 95 < 4.5, where a constraint penalizes a negative time derivative dq 95 /dt); (d) ℓ i3 (t); the vertical stability constraint from CREATE-NL is shown in black lines, the applied f margin factor is indicated in red text (note that two constraint curves are shown corresponding to the two different Ip traces, respectively the solid line of the initial trajectory and the dash-dotted line of the optimized trajectory).

Figure 12 .
Figure 12.A set of ramp-down simulations is presented to improve our understanding of the dynamics underlying the obtained optimum trajectory (for the case with dIp/dt = −100 kA s −1 , late HL transition, hot L-mode in figure 11).The initial, non-optimized ramp-down is shown in black solid lines.The case with optimized Ip and κ evolution is presented in green dash-dotted lines.A simulation applying the initial Ip evolution and the optimized κ evolution and a simulation with the optimized Ip evolution and the initial κ evolution are shown in blue dotted lines and red dashed lines respectively.(a) Ip(t); (b) ℓ i3 (t); (c) elongation κ(t); (d) q 95 (t); (e) ℓ i3 versus −Ip.

Figure 13 .
Figure 13.Initial and optimized time traces are shown for Ip, κ, q 95 and ℓ i3 (for the case with dIp/dt = −100 kA s −1 , late HL transition, hot L-mode in figure11).The values of plasma current and elongation of the equilibria initially obtained from CREATE-NL are indicated by the blue diamonds.The red diamond symbols indicate the equilibria calculated during post-processing of the optimized results.These equilibria are calculated with values of q 95 , κ and kinetic profiles consistent with the optimal trajectory.For an equilibrium at t = 55 s (identified by the green boxes): the initial equilibrium boundary shape obtained from CREATE-NL is compared to the boundary shape with optimized κ.(a) Ip(t); (b) q 95 (t); (c) elongation κ(t); (d) ℓ i3 ; (e) plasma boundary shape.
Figure 13.Initial and optimized time traces are shown for Ip, κ, q 95 and ℓ i3 (for the case with dIp/dt = −100 kA s −1 , late HL transition, hot L-mode in figure11).The values of plasma current and elongation of the equilibria initially obtained from CREATE-NL are indicated by the blue diamonds.The red diamond symbols indicate the equilibria calculated during post-processing of the optimized results.These equilibria are calculated with values of q 95 , κ and kinetic profiles consistent with the optimal trajectory.For an equilibrium at t = 55 s (identified by the green boxes): the initial equilibrium boundary shape obtained from CREATE-NL is compared to the boundary shape with optimized κ.(a) Ip(t); (b) q 95 (t); (c) elongation κ(t); (d) ℓ i3 ; (e) plasma boundary shape.

Figure A1 .
Figure A1.A comparison is shown between Te profiles for the DEMO stationary state calculated in RAPTOR, for different assumptions regarding the way radiated power is taken into account when calculating the loss power.The three profiles have an H factor equal to unity H 98y,2 ∼ (P L ) 0.69 τ E scl W th /(P L )τ E = 1.For the blue profile (P L )τ E = (P L )τ E scl = P heat − Ẇth , corresponding to the standard implementation of the gradient-based model as applied in this paper.For the red dashed profile (P L )τ E = (P L )τ E scl = P heat − Ẇth − P rad core , as proposed in [57].For the green dash-dotted profile, (P L )τ E scl = P heat − Ẇth − P rad core , while (P L )τ E = P heat − Ẇth .

Table 1 .
Operating parameters for stationary (flat loop voltage profile U pl (ρ)) flat-top burning plasma DEMO operating point, as evaluated by RAPTOR.

Table 3 .
Final (maximum) value of the internal inductance ℓ i3 for the RAPTOR simulations shown in figure 8. dIp/dt (kA s −1 ) t final (s) ,ℓ i3 ) plot in figure 10.For each of the considered rampdown rates dI p /dt = −50 kA s −1 , −100 kA s −1 , −150 kA s −1 , −200 kA s −1 , the two assumptions for the H 98y,2 time evolution introduced in section 6 are considered: early HL transition to cold L-mode versus late HL transition to hot L-mode.The maximum values of ℓ i3 reached at the final time of the simulation (for I p = 5 MA) are shown in table p 1,n1 , p 2,1 and p 2,n2 .The remaining parameters are assembled in the optimization vector p = [p 1,2 ...p 1,n1−1 p 2,2 ...p 2,n2−1 ] T .The optimization variables in this vector p contain the values of the actuator time traces on a set of knot points (the points indicated with star symbols in figure11, excluding the initial and final time points) and are varied in each iteration of the optimization routine.
P L )τ E (P L )τ E sclH 98y,2 (no P rad core ) Subtracting P rad core only in the scaling law evaluation would lead to the dash-dotted green profile in figureA1, with improved confinement.•Nostationaryoperatingpointwasfoundfor the case where P rad core is only subtracted for the confinement time evaluation.This can be understood by considering the corresponding expression:Evaluating this expression with the T e profile obtained with the standard settings of the gradient-based model (no P rad core subtraction, blue profile in figureA1), leads to an H factor above one (confinement is degraded, the blue T e is too optimistic).However, by decreasing T e , both numerator and denominator decrease, respectively due to a reduced W th and an increased P rad core .Interestingly, this leads to the fact that, when reducing T e (by reducing µ Te ), no profile satisfying H 98y,2 = 1 is found.In this paper, the H factor is calculated without subtractingPrad core in (P L ) τE and (P L ) τ E scl .Not applying the P rad core corrections as proposed in [57] in the gradient-based transport model is equivalent to assuming a confinement enhancement of about 10%.Van Mulders  https://orcid.org/0000-0003-3184-3361O Sauter  https://orcid.org/0000-0002-0099-6675E Fable  https://orcid.org/0000-0001-5019-9685F Felici  https://orcid.org/0000-0001-7585-376XM Mattei  https://orcid.org/0000-0001-7951-6584F Palermo  https://orcid.org/0000-0002-7524-3248A A Teplukhina  https://orcid.org/0000-0002-9213-7594