Axisymmetric predictions for mitigated and vertically unstable disruptions in ITER with runaway electrons

Abstract

We present two-dimensional global simulations of mitigated and vertically unstable disruptions in ITER in the presence of runaway electrons (REs). An elongated plasma in free-boundary equilibrium is subjected to an artificial thermal quench (TQ) and current profile flattening, followed by one or more massive material injections and a RE avalanche. Scenarios of major disruptions as well as upward and downward vertical displacement events are considered. Results provide important insights into the effects of RE formation, post TQ current profile, injection quantities and timings, and impurity flushout on the overall evolution of disruption and the plasma vertical motion thereof. Interplay between the various effects offers scope for potentially beneficial RE mitigation scenarios.

1 Introduction

Tokamak plasmas are known to occasionally get into states, wherein small perturbations grow exponentially leading to a disruption or termination of the discharge (Waddell et al. 1979). Such a disruptive event is typically described by various phases that the plasma undergoes. Nonlinear interaction of large-scale magneto-fluid instabilities cause a stochastization of the magnetic field structure, effectively causing most of the plasma thermal energy to be lost on a short time scale (${\sim }1$ ms in many existing devices). Such a thermal quench (TQ) phase increases the plasma electrical resistivity by orders of magnitude, and is therefore followed by a current quench (CQ) phase wherein the plasma current decays at a time scale of resistive diffusion. During the CQ phase, depending on the plasma conditions, it is possible that a seed of high energy relativistic electrons in the plasma grows exponentially so as to form a beam that carries nearly the entire plasma current by the end of CQ (Breizman et al. 2019). In fusion relevant tokamak devices, the current carried by such runaway electron (RE) beams can be a significant fraction of the pre-disruption current. Moreover, the preferred choice of vertically elongated plasmas in today's devices makes the plasma conducive to a vertical instability or a vertical motion of the plasma towards the first wall, often referred to as a vertical displacement event (VDE). A VDE can either precede and be the primary driver of a TQ (a hot VDE), or can follow the TQ in which case it is referred to as a cold VDE (Hender et al. 2007; Artola et al. 2024). While disruptions (in whichever form they occur) are not a serious issue in small and medium size tokamak devices, their consequences can be potentially severe in fusion-grade tokamaks such as the International Thermonuclear Experimental Reactor (ITER). They can manifest as one or more of potentially damaging thermal loads on the divertor, enormous electromagnetic forces on the surrounding conducting structures and RE induced thermal loads that can melt the first wall and beyond. It therefore becomes imperative to try to avoid or mitigate the effects of disruptions.

The aim of disruption mitigation is to simultaneously address all the detrimental consequences of disruptions: the thermal loads, electromagnetic loads and RE induced damage (Eidietis 2021). Towards this end, over the years several mitigation strategies have been proposed and evaluated for ITER. Currently, the ITER disruption mitigation system (DMS) is based on the injection of massive quantities of neon and/or deuterium into the plasma in the form of shattered pellets, in one or more stages (Lehnen et al. 2023). Depending on the specifics of any given disruption, there can be a few different schemes of executing the injections. Nevertheless, all the schemes essentially take advantage of the following general effects of neon and deuterium. Neon aids in effectively radiating the thermal energy, thereby reduces divertor thermal load and also avoids very long CQ times. In the context of RE beam mitigation, while impurities can cause a faster decay of the beam, they also aid in significantly increasing the avalanche growth (an unintentional effect in spite of the increase in critical electric field) via the additional bound electrons that act as additional targets for the avalanche process. On the other hand, deuterium during first injection can aid in plasma pre-dilution and reduced RE hot-tail seed (Nardon et al. 2020), and during a second injection onto a fully/partially formed RE beam, it can aid in neutralizing the ions (flushout) and cause a benign loss of REs (Reux et al. 2021; Hollmann et al. 2023). In addition to the above effects, inherent vertical instability of an elongated plasma (as is the case with ITER) leads to a vertical motion of the plasma column while it is going through massive material injections and RE formation. It is therefore evident that any reliable estimate or assessment of a disruption must self-consistently include the effects of RE formation, vertical motion and material injection at the very least. There have been attempts earlier towards such predictions mainly via the use of zero-dimensional (0-D) models (Kiramov & Breizman 2018; Martın-Solıs et al. 2022) using a force-free assumption along with a perfectly conducting vacuum vessel. The two-dimensional (2-D) equilibrium code DINA has also been employed to make similar predictions (Lukash & Khayrutdinov 2016), however, without considering the effects of partially ionized impurities with regards to RE source and effective electric field computations. Not considering spatial effects and the effects of partially ionized impurities on RE sources in a consistent manner can affect the predictions significantly. Furthermore, the aforesaid interplay can affect major disruptions (MDs), upwards VDE and downward VDEs very differently, and this can have important implications. This is the main motivation of the present work.

In this work we present free-boundary 2-D/axisymmetric simulations of an enlongated ITER-like plasma that undergoes a disruption, using the fusion magnetohydrodynamics (MHD) code JOREK (Czarny & Huysmans 2008; Hoelzl et al. 2021). The plasma is first subjected to an artificial TQ and current profile flattening, followed by a first injection of a mixture of neon and deuterium, alongside an introduction of an RE seed. Subsequent to this, the plasma goes through CQ and RE avalanche growth, vertical motion and a neon second injection and/or flushout. The interplay of the various effects are studied considering different disruption types (MD, upward and downward VDEs), different injection quantities and timings, and different types of current profile flattening. The effects of conducting structures surrounding the plasma are accounted for via the JOREK-STARWALL coupling (Merkel & Sempf 2006; Hoelzl et al. 2012; Artola 2018), while the REs are treated using an RE fluid model (Bandaru et al. 2024a) that is electromagnetically coupled to the background plasma. Effects of partially ionized impurities on the RE avalanche are taken into account through the models of Hesslow et al. (2018, 2019).

The paper is organized as follows. In § 2 we briefly present the physical model used for the simulations. In § 3 the details of the simulation set-up and the various physical properties used in the simulations are presented. This is followed by a detailed analysis of the results in § 4, followed by a summary of the key conclusions and outlook in § 5.

2 Physical model

As mentioned earlier, for the present simulations, we use the fusion MHD code JOREK. In particular, a reduced MHD version of the code is used that additionally accommodates impurities as well as REs. Both the deuterium ion mass density and total impurity mass density are evolved or solved for, among others, via time-dependent governing equations. However, individual charge states of the impurities are not evolved separately. Instead, the impurity charge state distribution is obtained using the coronal equilibrium model (Mosher 1974), while the code in general also has more advanced impurity models with higher computational costs (Hu et al. 2021). Material injection or flushout when needed is executed via appropriate density source or sink terms in the governing equations, and are activated temporarily for a short duration in order to raise or lower the densities to the required levels. The REs are treated as an extra fluid species separate from the background plasma and impurities. The RE fluid couples electromagnetically to the background plasma, via the modification that arises in the Ohm's law due to the REs not being subjected to resistive decay. Furthermore, only the density of the RE fluid ($n_r$) is evolved, without considering the dynamics of RE pitch angle and momentum. The RE density evolution includes an RE avalanche volumetric source along with parallel transport via a large parallel diffusivity and the $\boldsymbol {E}\times \boldsymbol {B}$ drift. Parallel diffusivity is used to mimic the effect of RE advection at the speed of light in a computationally cost effective manner (see Bandaru et al. 2019, 2021 for details). The RE avalanche source includes the effects of partially ionized impurities (i.e. the enhanced critical electric field as well as the bound electron targets for avalanche) through the models developed in Hesslow et al. (2018, 2019). However, the effect of impact ionization by REs on the charge state distribution of the impurities is neglected.

The electromagnetic coupling of all the passive electrically conducting structures as well as the current carrying coils are accounted for via the JOREK-STARWALL coupling. The coupling involves using non-local electromagnetic boundary conditions for the plasma using Green's function, without physically extending the domain or grid into the vacuum region. Furthermore, for the simulations presented in this paper, the diamagnetic drift as well as the field parallel velocity of the plasma fluid are not considered to reduce computational costs for the long time scales simulated in the present work. In other words, the background plasma fluid velocity is assumed to comprise only the $\boldsymbol {E}\times \boldsymbol {B}$ drift. Comprehensive details of the reduced MHD model as well as the numerical schemes involved in JOREK are described in Czarny & Huysmans (2008) and Hoelzl et al. (2021), while the specifics of the JOREK-STARWALL coupling can be found in Merkel & Sempf (2006), Hoelzl et al. (2012) and Artola (2018). The RE fluid model used hereby (along with model benchmark studies comparing with GO and DINA) have been described in detail in Bandaru et al. (2024a).

3 Overview of the simulation set-up

For the set-up used in this study, an elongated ITER plasma (pure deuterium) is considered that is in free-boundary equilibrium, with a toroidal plasma current $I_p=14.6$ MA, a uniform electron density of $n_e=10^{20}\ \mathrm {m}^{-3}$, an on-axis core electron temperature $T_{e,\mathrm {core}}\approx 15$ keV, central safety factor $q_0 \approx 1$ and an edge safety factor $q_\mathrm {edge}\approx 3.7$. The conducting structures (around the plasma) considered are the vacuum vessel (both the inner and outer shells), the poloidal field (PF) coils (PF1 to PF6), internal vertical stability coils (VS3), the central solenoid (CS1 to CS6), the outer triangular support (OTS) and the divertor inboard rail (DIR) (Artola et al. 2022). Among these, the vacuum vessel, OTS and DIR are passive structures, the most important in the present context of vertical stability in ITER being the vacuum vessel with an electrical resistivity $\eta _\mathrm {wall}=0.8 \times 10^{-6} \,\Omega \mathrm {m}$ and a thickness of $6$cm (each shell). The equilibrium plasma profiles and magnetic flux contours are shown respectively in figures 1(a) and 1(b), while the configuration of the conducting structures is shown in figure 1(c).

Figure 1. (a) Initial equilibrium plasma profiles of electron temperature $T_e$, toroidal current density ($J$), the safety factor ($q$) as a function of normalized poloidal magnetic flux $\psi _N$; (b) equilibrium magnetic flux contours; red: last closed flux-surface, blue: first wall, green: simulation domain used in JOREK; (c) Schematic showing the cross-section of the plasma along with the active and passive conducting structures of ITER used in the simulations. Note that the thick coils are discretized by several bands of thin structures in STARWALL. The acronyms used are defined as: VV for inner and outer shells of the vacuum vessel, PF for poloidal field coils, CS for central solenoid, D for divertor inboard rail and O for outer triangular support.

For the initial $1.5\ {\rm ms}$ in the simulation, the system is evolved from the initial static-equilibrium state resulting from the free-boundary Grad–Shafranov solver in JOREK to a state (that largely resembles the initial state) including the $\boldsymbol {E}\times \boldsymbol {B}$ flows that arise naturally.

At this time point, the various simulations can be seen as bifurcating into MDs and hot VDEs. For the MDs, an artificial TQ and current flattening is initiated immediately. For the hot VDEs, first a current perturbation is applied to the internal vertical stability coils (VS3) over a duration of $0.7\ {\rm ms}$. This causes a loss of vertical stability and a subsequent vertical motion either upwards or downward depending on the sign of the imposed perturbation. For these hot-VDE cases, the artificial TQ and current flattening is initiated $50\ {\rm ms}$ after the plasma axis has drifted by $\Delta Z_{a}=0.2\ {\rm m}$ from its initial equilibrium state. Such a timing corresponds to a detection of a hot VDE by the ITER control system within a vertical plasma movement of $0.2\ {\rm m}$ and a subsequent $50\ {\rm ms}$ lag for the arrival of impurities injected by the ITER DMS.

As mentioned earlier, at this time point, corresponding to $t=1.5\ {\rm ms}$ for MD and $t=t_{(Z_a=0.2m)}+50\ {\rm ms}$ for hot VDEs, the plasma is artificially cooled via a large perpendicular thermal diffusivity over a duration of $0.5\ {\rm ms}$ so as to attain an on-axis electron temperature $T_{e,\mathrm {core}} \approx 20\ {\rm eV}$. During this phase, the current density profile is simultaneously flattened artificially through the use of a large electrical hyper-resistivity (different levels of current profile flattening were explored). At the end of this phase, the artificial cooling is switched off, ohmic heating is switched on and a mixture of neon ($5\times 10^{21}$ or $1.5\times 10^{23}$ atoms) plus deuterium ($2\times 10^{23}$ atoms) as well as an RE seed population ($0.1$ A) is introduced within the whole plasma computational domain (including the halo or open field line region) with a spatially uniform density. This is intended to mimic an idealized massive material injection in the context of the ITER DMS. The injected neon and deuterium corresponds to a $10\,\%$ assimilation of a $28\ \mathrm {mm}$ ITER DMS pellet and a density rise of $\Delta n_\mathrm {Ne}=4.51\times 10^{18} \ \mathrm {m}^{-3}$ or $\Delta n_\mathrm {Ne}=1.35\times 10^{20} \ \mathrm {m}^{-3}$ and $\Delta n_{D}=1.81\times 10^{20}\ \mathrm {m}^{-3}$, respectively. Note that full simulations (until RE termination) of MD (up) and hot VDEs (up and down) are performed both with $\Delta n_\mathrm {Ne}=4.51\times 10^{18}\ \mathrm {m}^{-3}$ and a much larger quantity $\Delta n_\mathrm {Ne}=1.35\times 10^{20}\ \mathrm {m}^{-3}$ to explore the generation of RE beams at larger RE current. After the injection, due to impurity radiation, the plasma tends to further cool down to lower temperatures competing with ohmic heating. At the same time, due to the increased plasma resistivity at low temperatures, the plasma current decays (i.e. CQ), along with the conversion of thermal current to RE current via the avalanche mechanism. During and after the TQ, the plasma becomes further vertically unstable and moves upwards/downwards. This leads to a continuous shrinking of the plasma column due to scraping-off at the wall, and the eventual loss of the whole plasma.

Simulations are also performed with either a second injection of neon or a second deuterium injection via neon flushout for the MD and upward hot-VDE cases. This is executed either midway during the RE avalanche (when the RE current is roughly half the peak plateau RE current) or once the RE current reaches its plateau phase. The neon flushout is modelled via a density sink term (negative source). The list of simulations and their main differences are summarized in table 1. The runs with the larger quantity of neon are numbered with the suffix ‘L’ as shown in table 1. The downward hot-VDE case with $\Delta n_\mathrm {Ne}=4.51\times 10^{18} \ \mathrm {m}^{-3}$ does not produce any noticeable RE beam and so it does not make sense to perform a second neon injection or neon flushout for that case.

Table 1. List of simulations including MDs, up and downward hot VDEs. The third column refers to the neon density rise in the first injection during which the deuterium density rise is $\Delta n_{D}=1.81\times 10^{20}\ \mathrm {m}^3$. The fourth column refers to the change in the plasma internal inductance $\Delta l_i$ during the TQ phase that indicates the extent of current profile flattening. Terminology: Ne-2nd-P (second injection of neon in the early RE plateau phase); Ne-2nd-M (second injection of neon when about half the expected RE current is reached); D-2nd+Ne-Flush-P (second injection of deuterium along with complete flushout of neon starting in the early RE plateau phase); D-2nd+Ne-Flush-M (second injection of deuterium along with complete flushout of neon starting when about half the expected RE current is reached). Run names with the suffix ‘L’ indicate a 30 times larger neon first injection quantity than those without the suffix.

3.1 Further details of the simulations

The artificial or pre-first-injection cooling of the plasma (e.g. from $t=1.5\ {\rm ms}$ to $t=2\ {\rm ms}$ in MD simulations) is achieved by the use of a large perpendicular thermal diffusivity. During this phase of the simulation, the parallel thermal diffusivity is increased simultaneously, such that there is no accumulation of thermal energy lost from the plasma core in the scrape-off layer (SOL) and the resistivity consequently remains high there. This ensures that the SOL is relatively cooler and the halo currents are minimized. The choice to minimize halo currents was made in order to maximize RE generation, since the current spread into the halo region could decrease the available electric field in the confined plasma. The plasma electrical resistivity ($\eta$) used is a function of temperature and $Z_\mathrm {eff}$, and is given by $\eta _\mathrm {sp} = \eta _{\mathrm {sp},c} \times Z_\mathrm {eff}$, where $\eta _{\mathrm {sp},c}$ is the classical Spitzer resistivity in the absence of impurities. In this respect, the dependence of the Coulomb logarithm on temperature and density is not considered and instead a constant value ($\ln \varLambda =20$) has been used that corresponds to a hot plasma at $T_e=15\ {\rm KeV}$ and $n_e=2\times 10^{20}\ \mathrm {m}^{-3}$. The background plasma fluid viscosity ($\mu =5.2\times 10^{-7}\ {\rm kg}\ ({\rm m}\ {\rm s}^{-1})$) and parallel thermal diffusivity ($\kappa _\parallel = 1.54 \times 10^{29}\ \mathrm {m}^{-1}\ \mathrm {s}^{-1}$) are chosen to be temperature independent, and there are no volumetric sources of thermal energy besides ohmic heating and radiation losses. The SOL diffusivities are not treated differently. Standard fixed boundary conditions are used for the temperature, with the far SOL $T_e\approx 2\ {\rm eV}$. For the electric potential, we used the Dirichlet condition. Sheath boundary conditions might be more realistic but is beyond the scope of the present work.

A spatially uniform RE seed current of $0.1$ A has been used throughout and intends to represent the primary seed of REs that would have been generated during the TQ via the mechanisms of hot-tail, Compton scattering, tritium decay and Dreicer (Breizman et al. 2019). While a significant variation in the magnitude of the RE seed is possible in reality, the plateau RE current has been observed to have a very weak dependence on the RE seed (Bandaru et al. 2024a), justifying the simpler choice made here. Likewise, the profile of the RE seed in reality depends on the precise details of the TQ phase involving flux-surface breakup and reformation along with current profile flattening. Precise RE seed profiles can only be obtained through comprehensive three-dimensional (3-D) MHD simulations of the TQ phase, which is beyond the scope of the present work. To reduce computational costs for the long time scales simulated here, RE parallel transport is modelled via a large parallel diffusivity (to mimic the fast parallel advection). A value of $D_{\mathrm {RE},\|}=1.54\times 10^{6}\ \mathrm {m}^2\ \mathrm {s}^{-1}$ has been used for the results presented in this work, along with a small perpendicular diffusivity of $D_{\mathrm {RE},\perp }=4.6\times 10^{-2}\ \mathrm {m}^2\ \mathrm {s}^{-1}$ for numerical stability. The choice for this value has been arrived through a sensitivity study that showed that much higher values would cause artificial perpendicular diffusion of REs in the later stages of the simulation. The perpendicular transport of REs occurs through the $\boldsymbol {E}\times \boldsymbol {B}$ advection of the RE fluid.

As mentioned earlier, for the sake of simplicity, deuterium neutrals are not considered. This is justified by the fact that for the temperatures involved in this simulation, the neutral population is expected to be small and will have an insignificant effect on the final results. This has indeed been confirmed by comparing with a test simulation involving deuterium neutrals. Also, we do not consider the parallel velocity of the background plasma in these simulations. Several test simulations showed that use of the customary fixed boundary condition leads to large density gradients near the domain boundary leading to a significant loss of particles over the long time scale of the simulation. To avoid this, natural or Neumann conditions have been used for the main-ion (deuterium) density and impurity (neon) density evolution equations, which ensures strict particle conservation in the domain. This is due to the choice of vanishing of flows normal to the domain boundary. Along with this, we use large values of particle diffusivity to maintain a nearly uniform density distribution. Values of $D_{\perp }=15.4\ \mathrm {m}^2\ \mathrm {s}^{-1}$ and $D_{\parallel }=1.54\times 10^{5}\ \mathrm {m}^2\ \mathrm {s}^{-1}$ respectively have been used for both the main ions and impurities.

A polar grid with radial and poloidal grid resolution $N_r \times N_\theta = 101\times 128$ bi-cubic Bezier finite elements is used with poloidal grid clustering in the region covered by the plasma during its vertical trajectory. We now turn to the discussion and analysis of the results obtained from the simulations.

4 Results and discussion

4.1 General characteristics and the effect of RE beam formation

Firstly we look into the general physical behaviour in our simulations, wherein the focus is on how RE formation modifies the overall evolution of a disruption. As mentioned earlier, the artificial TQ executed over a time of ${\approx }0.5\ {\rm ms}$ causes the on-axis temperature to drop to $T_{e,\mathrm {core}} \approx 20\ {\rm eV}$. After the first injection of deuterium and neon, impurity radiation dominates over ohmic heating causing a further cooling of the plasma to settle at roughly $T_{e,\mathrm {core}} \approx 10\ {\rm eV}$. Subsequently RE formation causes cooling of the plasma to continue to further lower temperatures due to the decrease in thermal plasma current and associated ohmic heating, such that the effective ion charge is around unity even in the presence of impurities. Figure 2 compares the electron temperature profile after the full RE beam formation for a MD case (run 1) to that when the RE avalanche is switched off (run 0). Additional cooling due to RE beam formation is evident. Moreover, the formation of RE beam also results in an increase in plasma internal inductance $l_i$ due to the general tendency of a peaked RE current profile formation. This can in principle affect the kink stability of the current beam, e.g. due to a relatively lower central safety factor.

Figure 2. Electron temperature profiles before and after RE formation for run 1 compared with the corresponding case wherein REs are switched off (run 0).

Figure 3 shows the subsequent effects of RE formation as compared with a situation when REs are switched off. It can be observed that the primary effect of RE beam formation is a decrease in the current decay rate post the beam formation. This is due to the fact that REs are not subjected to a decay on the resistive time scale. This also results in a relative slowing down of the plasma vertical motion towards the wall, since the current decay rate is a significant driving cause for the plasma vertical motion in ITER due to the relatively long wall time ${\sim }500\ {\rm ms}$. The slower plasma current decay rate with REs results in the corresponding slower increase in the toroidal wall current $I_w$ (see figure 3d) and also a relatively more gradual rise in the toroidal halo currents (see figure 3b). However, the total plasma current is roughly the same in both cases at the time when the current beam collapses into the wall. This can be clearly seen in figure 4, wherein the currents at any given vertical location of the plasma axis are plotted. Therefore, RE beam formation makes the plasma take longer to reach the wall (an additional ${\approx }23\ {\rm ms}$ in this specific case), but with roughly the same total plasma current. Qualitatively this aspect agrees well with the previous predictions of (Kiramov & Breizman 2018; Martın-Solıs et al. 2022). As mentioned earlier, the temperature of the SOL is maintained at ${\approx }2\ {\rm eV}$ via a Dirichlet boundary condition and a large parallel thermal diffusivity. The resulting higher plasma resistivity in the SOL relative to the plasma core ensures that the toroidal halo currents do not grow to large values relative to the total plasma current (except when it nears a complete scrape-off). This is a simplification made in the present study and a common feature in the simulations presented in this work. While a more comprehensive treatment of the SOL would definitely be more predictive, it is beyond the scope of the present work.

Figure 3. Comparison of the evolution of 0-D parameters during the MD cases with REs (run 1, black) and with REs switched off (run 0, orange). Nomenclature: $I_p$ is the total plasma current, $I_{\mathrm {pol},\mathrm {halo}}$ and $I_{\mathrm {tor},\mathrm {halo}}$ are the poloidal and toroidal halo currents, respectively, $I_\mathrm {RE}$ is the RE current, $q_{95}$ is the safety factor at $\psi _N=0.95$, $Z_a$ is the vertical position of the plasma magnetic axis, $l_i$ the plasma internal inductance, $I_w$ the wall current and $F_z$ the total vertical force on the vacuum vessel. Densities plotted in the panel (a) are average values over the whole simulation plasma domain (including the halo region) and so they remain constant even during scraping-off due to the expansion of the halo region.

Figure 4. Total plasma current ($I_p$) and RE current ($I_\mathrm {RE}$) versus the vertical position of the plasma magnetic axis ($Z_{a}$) for the MD cases with REs (run 1) and without REs (run 0).

Figure 5 shows the profiles of current density for run 1 at different time points until the RE plateau formation. The J profile after avalanche is a result of both the avalanche as well as the minor radial diffusion of the parallel electric field during that phase. One can observe that the J profile after full RE beam formation (green curve) is centrally peaked. Note that the magnitude of $J$ increases due to the plasma shrinking already during this phase. In an axisymmetric situation, one relevant/important force on the vacuum vessel is the vertical force $F_z$. However, it is important to note that the force on the vacuum vessel peaks only after the change in magnetic field diffuses into the vessel material, which occurs over $\tau _w = L/R \sim 500\ {\rm ms}$ time scale (Pustovitov, Rubinacci & Villone 2017; Clauser, Jardin & Ferraro 2019; Artola et al. 2022). The simulations performed here, although already on a relatively long time scale, only span until the  final collapse of the current beam that is ${\sim }100\unicode{x2013}150\ {\rm ms}$. This does not allow for a peaking of the wall force, and hence, the force values seen should not be interpreted as the maximum forces predicted.

Figure 5. Profiles of plasma current density for run 1 (MD) at various time points from the end of TQ till RE plateau formation. Bold black: immediately after artificial TQ and current flattening; dotted black: slightly after the end of the first injection of neon and deuterium; magenta: just before any significant RE beam formation; green: after full RE beam formation.

4.2 Impact of the post TQ current density profile

It is well established that the aspect of magnetic helicity conservation during the fast TQ phase of a disruption tends to flatten the current density profile relative to the pre-TQ profile (Biskamp 1993). In this regard, we explore the effect of the immediate post-TQ current density profile on further evolution of the disruption. What might be the most realistic current flattening scenario for a real disruption would likely depend on the details of MHD activity and magnetic stochastization before and during the TQ. For lack of known constraints on the post-TQ current profile, a few different current profile flattening scenarios during the artificial TQ phase have been chosen in the context of MDs. They include the cases of significant flattening (run 1), full flattening (run 7) and full flattening but only within the region enclosed by the $q=2$ flux surface (run 8). The corresponding post-TQ current density profiles are shown in figure 6(a), wherein they are characterized also by the change in internal inductance $\Delta l_i$ after the flattening.

Figure 6. (a) Profiles of plasma current density with different levels of flattening immediately after the artificial TQ ($t=2\ {\rm ms}$) as a function of normalized poloidal flux. (b) Relaxation of post TQ flattened profiles of plasma current density as a function of normalized poloidal flux. Black lines: run 1; magenta lines: run 7; blue lines: run 8.

The differences in the overall evolution in these cases is shown in figures 7 and 8. We observe a clear effect on the well-known $I_p$ spike (temporary rise in plasma current) that is often seen during and immediately after the TQ of a plasma disruption. The extent of the $I_p$ spike is upto $0.9$ MA in the cases considered and is observed to be directly correlated to the level of current profile flattening that the plasma is subjected to during the TQ phase (strong flattening leading to a larger spike). Furthermore, it can be seen that except for a short time window of ${\sim }10\ {\rm ms}$ post TQ, there is negligible difference in the evolution of run 1 (black) and run 8 (blue). This can be attributed to the relaxation (via resistive diffusion) of the current density profiles after the flattening. Both the profiles develop large near-edge gradients in the current density immediately after flattening (figure 6a). However, after the artificial TQ and the subsequent first injection of neon, the plasma resistivity is rather large and this ensures that the near-edge gradients are smoothed out within ${\sim} 20\ {\rm ms}$. This is shown in figure 6(b), where in fact by $t\approx 22.5\ {\rm ms}$, the core J profiles of all the cases are very similar, and in particular the profiles of run 1 and run 8 match very closely.

Figure 7. Comparison of the evolution of 0-D parameters for the MD cases with different current profile flattening during the artificial TQ. Due to fast current profile relaxation post flattening, the difference in evolution between run 1 and run 8 is negligible, while run 7 evolves differently due to different wall current distribution.

Figure 8. Total plasma current and RE current versus the plasma axis vertical position for the MD cases with different current profile flattening during the artificial TQ.

In spite of this effect, interestingly, the evolution of the case with a fully flat J profile (run 7) shows significant differences over the course of the simulation. This is due to a different distribution of current density on the plasma, which leads to a different wall current distribution when the plasma current decays. Such current density distributions delay the vertical movement (it even moves downwards during the initial phase after reversing its direction to upwards). One can see run 7 to be a case of being near the threshold between upward and downward motion. In summary, while the precise post-TQ J profile can have an effect on the further evolution of the disruption, fast radial resistive diffusion of current tends to suppress the differences to a great extent. We now turn to the effects of a second injection of massive material into the plasma.

4.3 Neon second injection and neon flushout

4.3.1 Neon secondary injection

A second injection of massive quantities of neon through a shattered pellet injection (SPI) onto an existing RE beam is a strategy to ensure a faster decay of the RE beam current. In our simulations, for the neon second injection cases, a total of $1.5\times 10^{24}$ neon atoms are introduced uniformly into the computational domain (assuming a $100\,\%$ assimilation and corresponds to a neon density rise $\Delta n_\mathrm {Ne}=1.35\times 10^{21} \ \mathrm {m}^{-3}$). This is done via an impurity particle source activated over a relatively short duration. The injected neon quantity amounts to a factor of $300$ times the existing neon atoms from the first injection, and about $5.4$ times the total existing deuterium atoms in the plasma (see figure 9a). The relatively large quantity of neon for second injection has been chosen intentionally to assess whether any beneficial effects could be expected at all. For the MD cases (branched from run 1), a second neon injection is made both at the beginning of the RE plateau phase (run 2) and at midway during RE avalanche i.e. when the RE current is half the expected peak value (run 4). For the upward-VDE case (branched from run 9), a second neon injection is made only at the beginning of the RE plateau phase (run 10). It must be noted that before a neon second injection is made, the plasma is already cold (${<}10\ {\rm eV}$) and so only the lower charge states of neon exist, meaning that $Z_\mathrm {eff}\approx 1$. So an additional neon injection does not alter $Z_\mathrm {eff}$ much and, therefore, the plasma electrical resistivity $\eta$ does not change significantly. It must be noted that neon injection executed here as a spatially uniform source over a short time scale is rather idealized. Clearly, material assimilation during an SPI in ITER will be via neon pellet penetration, ablation and transport that is spatially strongly non-uniform and varying with time. Such a complex process cannot be treated in a comprehensive manner in the present context of 2-D simulations, which is a limitation of the present study. Nevertheless, the uniform injection hereby aims to simplify the assimilation process in order to obtain an understanding of the global picture of plasma evolution after such an injection.

Figure 9. Comparison of the evolution of 0-D parameters for the MD cases with a neon second injection and neon flushout with a deuterium second injection both starting at the plateau phase and midway during the RE avalanche.

Comparison of the effect of a second neon injection in MD cases (runs 1, 2 and 4) is shown in figures 9 and 11(a). For the case with plateau injection (run 2), introduction of massive quantities of neon causes a significantly faster decay of the RE beam current, as can be seen from figure 9, due to the increase in the effective critical electric field. This correspondingly decreases the time scale of the vertical motion and so the plasma dumps into the wall much faster too, as compared with the case without a second injection. In this case, the net effect is that the RE current just before the final loss of the confined plasma region is lower, i.e. ${\sim} 1$ MA vs $1.5$ MA, for the case without a second injection, as can be seen from figure 11(a). The effect is similar for the case of the plateau neon second injection onto the upward VDE (run 10). However, in the upward-VDE case, the RE current before final collapse is only marginally lower than the corresponding case without a second neon injection (see figure 11b). This can be attributed to a slightly lower stabilization factor $f$ (Leuer 1989) in the case of an upward VDE, which determines the growth rate of the vertical motion.

Figure 10. Comparison of the evolution of 0-D parameters for the upward-VDE cases with a neon second injection and neon flushout with a deuterium second injection.

Figure 11. Total plasma current and RE current versus the plasma axis vertical position for (a) MD cases with a neon second injection and neon flushout with a deuterium second injection starting both at the plateau phase and midway during the RE avalanche. (b) Upward-VDE cases with a neon second injection and neon flushout with a deuterium second injection starting at the RE plateau phase.

When neon is introduced halfway during the avalanche growth (run 4), the immediate effect is quite different in the sense that the RE avalanche growth rate increases dramatically. This is due to the fact that there is enough electric field left and the number of electrons available for avalanche (which also includes bound electrons) increases drastically. That is, the effect of additional bound electrons from the added neon that aid as targets dominates the effect of the corresponding higher effective critical field $E_{c,\mathrm {eff}}$. The peak RE current obtained is however the same. It must be noted that the propensity to avalanche (due to a larger neon quantity and thereby more bound electrons) does not always imply that the final RE beam current should be higher. The parallel electric field falls due to RE avalanche and resistive current decay, and so the avalanche growth stalls once there is not enough electric field left. Further from the RE plateau downstream, the behaviour is similar to the other second neon injection cases, i.e. faster RE current decay and vertical motion.

In summary, a neon second injection in general causes faster RE beam decay and correspondingly faster vertical motion to the wall. Earlier injection (during avalanche) enhances the avalanche rate but the subsequent evolution is similar to the plateau injection.

4.3.2 Neon flushout induced by a deuterium secondary injection

In recent years, impurity flushout through a massive secondary injection of deuterium has been found to show some promise in obtaining benign RE beam terminations (equal to no wall damage due to REs) in experiments at DIII-D and JET (Paz-Soldan et al. 2021; Reux et al. 2021). It appears that a benign termination after neon flushout is caused not by a single or unique mechanism. Some known pathways to such a benign termination post neon flushout are hollow current density profile formation that are inherently MHD unstable (Reux et al. 2021) and a drop in $q_{95}<2$ but with a monotonous J profile (Paz-Soldan et al. 2021) etc. All these pathways are presumed to result in a violent MHD activity that distributes the REs relatively broadly over the wall. Importantly, the experimental observation seems to be rather robust. In experiments, after a deuterium second injection, impurity flushout from the plasma occurs due to neutralization and subsequent outward transport of the neutral impurity atoms. In our simulations, we model neon flushout in a simplified way via an artificial neon sink (negative source) to ensure an exponential loss of impurities on a $10\ {\rm ms}$ time scale (only about $1\,\%$ of impurities remain in the plasma after $10\ {\rm ms}$). The time scale of impurity flushout in DIII-D and JET was observed to be ${\sim }10\ {\rm ms}$. While there is no established consensus on what would be the expected neon flushout time scale in ITER, it can be expected to be higher due to the relatively larger spatial scales. Nevertheless, our choice of the time scale $10\ {\rm ms}$ is driven by the fact that the total time to RE beam loss itself is at most $50\ {\rm ms}$, and so a flushout on a much slower time scale would not have a significant effect. A total of $2.5\times 10^{24}$ deuterium atoms are injected uniformly in the domain, which amounts to about $9$ times the existing deuterium atoms in the plasma.

In the present simulations, the effect of a deuterium second injection plus a neon flushout is shown in figures 9–11, for both the MD cases and upward hot VDE (runs 5, 6 and 11). In the MD cases (runs 5 and 6) the impact is not as significant as a neon second injection and the result is a marginally faster RE current decay and corresponding vertical motion (see figure 9). In spite of a loss of neon impurities, we see a faster RE current decay. This is due to the fact that the pre-existing quantity of neon is relatively not large and the RE current decay is mainly dominated by the scrape-off rather than due to impurities. Hence, the loss of impurities is offset by the massive quantity of deuterium atoms injected. Thus, the net effect is a slightly faster decay of RE current. In a situation when the time scale of neon flushout is much larger, the effect of faster decay of the RE current could start to occur with a significant delay and can possibly even not show up before the RE beam loss. It must be noted that in many experiments a slow rise in RE current is observed post neon flushout; e.g. shot:95135 of JET (Reux et al. 2021)). Such an RE current rise can be attributed to the fact that in many such RE experiments, there is no scrape-off effect since the plasma is vertically stable (nearly circular limiter configuration).

In summary, in comparison with a neon second injection, neon flushout has a qualitatively similar but relatively weaker effect on the overall VDE. However, unlike in many existing experiments with circular plasmas, neon flushout in elongated plasmas can in fact lead to a current decay rather than a current rise, mediated by the scrape-off effect.

4.3.3 Estimates for RE energy loss after a neon second injection and flushout

In the context of massive material injections, it is important to assess how much net energy will have been extracted by REs that eventually gets dumped onto the first wall. While REs gain energy via acceleration by the parallel electric field they lose energy via the channels of synchrotron radiation, Bremsstrahlung and collisional drag. The synchrotron loss (per particle) is strongly dependent on the energy and pitch angle of the REs and scales (for small pitch angles) as $\propto B^2 \theta ^2 \gamma ^2$, where $B$ is the magnetic field, $\theta$ the RE pitch angle and $\gamma$ is the relativistic factor. On the other hand, the Bremsstrahlung and collisional drag losses per particle scale as $\propto n_e \gamma (Z_\mathrm {eff} + 1)$ and $\propto n_e (Z_\mathrm {eff} + 1 + \gamma ) \gamma ^{-1}$, respectively, again assuming a small pitch angle and high energy of the REs. Clearly, as mentioned earlier, both the pitch angle as well as the energy of the REs are not tracked in the RE fluid model used in the present simulations (note that future work including REs in a particle in cell (PIC) treatment in JOREK will improve the modelling fidelity). Nevertheless, it is still possible to make indirect rough estimates of the RE energy losses via the effective critical electric field $E_{c,\mathrm {eff}}$ that is computed within the RE fluid model while evaluating the avalanche source. That an estimate can be made using $E_{c,\mathrm {eff}}$ would be evident from the fact that the analytical model for $E_{c,\mathrm {eff}}$ used in our RE fluid model as derived by Hesslow et al. (2018) already includes the effects of synchrotron radiation, Bremsstrahlung and collisional loss. However, it must be noted that due to the nature of the assumed distribution, the model of Hesslow et al. (2018) might not be accurate at electric fields that are much higher than the critical electric field. This ($E \gg E_{c,\mathrm {eff}}$) can happen in certain phases of the simulation (e.g. during avalanche), which can affect the accuracy of the estimates made hereby, in terms of underestimating the synchrotron losses during such a phase. Nevertheless, this would imply an overestimate of the undissipated RE energy values presented hereby and would only be on the conservative side.

We now proceed to make RE energy loss estimates as follows. The parallel momentum balance for REs can be written as

(4.1)\begin{equation} \frac{{\rm d}}{{\rm d}t} \left(\gamma m_{e0} v_{\|} n_r \right) ={-}en_r E_{\|} - F_\mathrm{drag}. \end{equation}

Multiplying the above equation with $v_{\|}$, using $F_\mathrm {drag}=-e n_r E_c^\mathrm {eff}$ and integrating over the computational volume, one can rewrite the above equation as

(4.2)\begin{equation} \int \frac{{\rm d}}{{\rm d}t} \left(\frac{1}{2}\gamma m_{e0} v_{\|}^2 n_r \right) \,{\rm d}V\, {\rm d}t = \int J_{RE,\|} \left(E_{\|}-E_c^\mathrm{eff}\right) \,{\rm d}V \,{\rm d}t, \end{equation}

So the total poloidal magnetic energy channelled to REs is $\int J_{RE,\|} E_{\|} dV dt$ (which is the total work done by the electric field on the REs), out of which $\int J_{RE,\|} (E_{\|}-E_c^\mathrm {eff})\,{\rm d}V\,{\rm d}t$ goes towards their kinetic energy and $\int J_{RE,\|} E_c^\mathrm {eff}\,{\rm d}V\,{\rm d}t$ is the energy needed to sustain the RE current. At the end, this total energy channeled to the REs gets lost by collisions, Bremsstrahlung, synchrotron and deposition on plasma facing components (PFCs) via scrape-off and final loss.

The total energy to REs in different injection scenarios is plotted in figure 12 for the MD cases. Expectedly, both the cumulative energy channeled to REs (by the parallel electric field) as well as the cumulative energy dissipated by them (via collisions, Bremsstrahlung and synchrotron) increase over time. However, it can be observed that in all the cases shown, nearly $70\unicode{x2013}100$ MJ of the total energy extracted by the REs is undissipated by the time of final termination (the difference between the bold lines and dashed lines in figure 12). This can be attributed to the fact that while a second injection causes an increase in RE energy dissipation, at the same time it also increases the total energy extracted by the REs from the poloidal magnetic field. Therefore, it turns out that the net effect of a second neon injection and neon flushout therefore is rather marginal, when it comes to the cumulative undissipated energy of the REs by the time of final termination. It must be noted that the aforementioned undissipated RE energy not only includes the energy remaining with the REs at final termination, but also the total RE energy lost via scrape-off of REs over the whole duration of the simulation. The present results are qualitatively different from the conclusions arrived by Martın-Solıs et al. (2022), wherein neon injection is seen to significantly decrease the undissipated RE energy (and the decrease is proportional to the amount of neon injected). However, our observation here is qualitatively similar to the conclusion arrived by earlier work using DINA (Lukash & Khayrutdinov 2016), though there are several differences between both the models (e.g. DINA's model did not include the effects of partially ionized impurities). Similar to the enhancement of total RE beam current formed after a neon first injection (due to additional target electrons), sustained avalanche post neon second injection (runs 2 and 4) or a deuterium second injection/neon flushout (runs 5 and 6) ensures a slightly larger energy extracted by the REs.

Figure 12. Cumulative poloidal magnetic energy that is channeled to REs $\int J_{RE,\parallel } E_{\parallel } \,{\rm d}V \,{\rm d}t$ (bold lines) and the part $\int J_{RE,\parallel } E_{c, \mathrm {eff}} \,{\rm d}V \,{\rm d}t$ that is dissipated by REs via collisions, synchrotron and Bremsstrahlung (dashed lines), in different injection scenarios. The difference between the bold and dashed curves indicates the instantaneous RE energy content plus the cumulative RE scrape-off energy loss until that point in time.

Figure 13. Normalized toroidal current density distribution just before the initiation of artificial TQ for the (b) MD (run 7), (c) up hot-VDE (run 9) and (d) down hot-VDE (run 12) cases.

Figure 14. Comparison of the evolution of 0-D parameters comparing the baseline cases (without second injections) of MD (run 7) with that of the hot-VDE up (run 9) and hot-VDE down (run 12). Note that the sign for the wall force for the downward VDE case (run 12) has been reversed for easy comparison.

Figure 15. Total plasma current and RE current versus the vertical distance of the plasma magnetic axis from the wall touching point, comparing the baseline cases (without second injections) of MD (run 7) with that of the hot-VDE up (run 9) and hot-VDE down (run 12).

Figure 16. Maximum RE beam current formed versus the first injection quantity of neon. The injected deuterium was fixed at $n_{D}=1.81\times 10^{20}\ \mathrm {m}^{-3}$.

Figure 17. Comparison of the evolution of 0-D parameters for the baseline MD cases with different first injection neon quantities (run 1 and run 1L).

Figure 18. Total plasma current and RE current versus the plasma axis vertical position (run 1 and run 1L).

4.4 Major disruptions vs hot VDEs

In this section we look into the main differences between the MD (run 7) and hot-VDE baseline cases (runs 9 and 12), all without any additional secondary injection of material. In the hot-VDE cases, effectively the artificial TQ gets initiated about ${\sim }75\ {\rm ms}$ later than that in the MD case. During this time, due to the initially applied current perturbation in the internal vertical stability coils (VS3), the plasma in the hot-VDE cases drifts significantly (and reduces in size) already before the TQ is initiated. This is evident from the vertical parts of the dashed blue and magenta curves at nearly constant $I_p$ in figure 15, wherein the axis distance from the wall is plotted versus the total and RE currents. The relative plasma positions just before the TQ is shown in figure 13 for all the three cases.

In the case of the upward hot VDE (run 9), the overall evolution of the disruption is in general qualitatively similar to the MD case (run 7). There are however certain differences, that arise due to the delay in TQ initiation in the hot VDE. A relatively larger $I_p$ spike is observed in run 9 as compared with run 7 due to the slightly higher flattening of the current profile during the artificial TQ. Furthermore, due to the relatively smaller cross-section by the time of first injection, the total assimilated neon within the last closed flux surface (LCFS) is smaller in this case (since the total number of injected atoms remains the same). As a consequence, the avalanche-magnifying effect of partially ionized impurities is reduced and, therefore, results in a relatively smaller plateau RE current of ${\sim }3.7$ MA vs ${\sim }5$ MA in the case of MD. However, the later formation of RE current in run 9 eventually leads to a similar RE current at the time of final collapse (${\sim }1$ MA), as is evident from figure 15.

The situation with the downward hot-VDE case (run 12) is however starkly different. Although the pre-TQ VDE growth rate is similar to that of the upward hot VDE (run 9), it gets much higher post the TQ. This can be seen from figures 14(d) and 15. This leads to a much faster shrinking of the plasma such that the effect of scraping-off dominates the RE avalanche, consequently leading to nearly no RE beam formation at all. In addition, in this situation, the plasma confined region vanishes at a significantly larger plasma current (see figure 15). The faster time scale of vertical motion (and plasma shrinking) even at larger $I_p$ also leads to higher halo currents as can be seen from figure 14. Faster post-TQ vertical motion in the downward VDE can be attributed to the relatively stronger gradient in the magnetic field below the lower X point in ITER, that in turn causes a stronger vertical force on the plasma during the downward VDE. Nevertheless, such a scrape-off dominated scenario could point to a potential route for the avoidance of RE beam formation albeit with the consequence of a higher total current at final collapse and larger wall forces. Prior work (Lukash & Khayrutdinov 2016) using DINA simulations of ITER downward hot VDEs showed a similar behaviour, wherein the vertical motion and scrape-off was relatively far stronger than in upward VDEs. While it can be difficult to control the direction of the plasma vertical motion in ITER, in scenarios wherein downward hot VDEs occur naturally, such an RE-free scenario can be potentially beneficial.

4.5 Effect of a first injection neon quantity

The effect of a first injection quantity of neon on the RE beam current formed in MD cases based on run 1 is shown in figure 16. Injection of lower quantities of neon (first injection) reduces the RE beam current formed. In fact, with a neon injection quantity that is $3\,\%$ of that in run 1, it has been observed that there is no noticeable RE formation even after $100\ {\rm ms}$ (run 1a). This is due to the relatively quick and significant ohmic reheating of the plasma in the absence of much neon, post which there is not enough electric field available for an efficient RE avalanche. The present observation is in line with previous predictions with partially ionized impurities (Vallhagen et al. 2020) and also with the usual absence of significant post-disruption RE beams in existing experiments without high-Z injections. However, with present assumptions in ITER, the requirement to radiate at least ${\sim }85\,\%$ of thermal energy and to maintain a CQ time $t_\mathrm {CQ}<150\ \mathrm {ms}$ would need a neon first injection density rise ${\sim} 10^{19}$ or higher (it is possible that the required rise can be slightly lower in case the assimilation of hydrogen from the first injection is higher). This emphasises that the risk from REs in ITER is strongly enhanced by high-Z injections meant to simultaneously mitigate high thermal loads on the divertor and high electromagnetic loads on the vacuum vessel. It must be noted however, that increasing the injected neon beyond a certain level does not translate to a commensurate rise in the RE beam current (the extent of electric field available becomes important and limits the obtained RE beam current).

Figure 19. Comparison of the evolution of 0-D parameters for the baseline downward-VDE cases with different first injection neon quantities (run 12 and run 12L). Note that the label $Z_a$ in this specific case (downward motion) represents $-Z_a$.

Figure 20. Total plasma current and RE current versus the plasma axis vertical position (run 12 and run 12L).

4.5.1 Effect of a neon first injection quantity on the overall evolution of MD and VDE cases

As mentioned earlier, full simulations for MD (up), upward VDE and downward VDE have been performed also with a neon first injection quantity that is $30$ times larger. These simulations are referred to with the suffix ‘L’ in the table 1. Here we compare the effect of the much larger neon injection on the evolution of the disruption.

The overall evolution is in many ways qualitatively similar even for the larger neon injection. However, there are differences in details that arise. Figure 17 highlights the important differences in evolution of the global parameters for the MD baseline cases run 1 and run 1L. As expected, a higher neon quantity leads to a lower temperature after first injection, and therefore, a faster CQ. This causes the $q_{95}$ to peak much earlier since this is a phase when current decays without the plasma shrinking in size significantly ($q_{95} \sim a^2/I_p$). Faster CQ also causes the threshold $I_p$ value (${\sim }10$ MA) to be reached earlier and, hence, the plasma vertical instability starts earlier in the case of a larger neon injection. Likewise, an increase in the number of bound electrons available for avalanche leads to a larger RE beam current ${\sim }9.5$ MA at a higher neon injection (run 1L) versus ${\sim }5.5$ MA for run 1. Significant plasma vertical motion and the associated halo current ${\sim }3\,\mathrm {MA}$ already by the time of RE avalanching in run 1 also contributes to this difference in final RE current. However, more important is the decay rate of the RE beam current. In the case of run 1L, when the full RE beam current is formed, it is not immediately subjected to a decay via scrape-off unlike in run 1. Due to this the plateau phase is relatively flat and so the plasma motion toward the wall is much slower. Therefore, unlike in run 1, in the case of larger neon injection, the decay of $q_{95}$ after peaking is a result of plasma shrinking at a nearly constant $I_p$. This can be seen clearly in figure 18, in the vertical leg of the $I_p$ curves. However, what is practically important eventually is the fact that with a larger neon injection, due to relatively prompt conversion via RE avalanche, $q_{95}<2$ occurs far into the RE plateau phase. Differences in the overall evolution described here is quite similar in the case of upward-VDE cases run 9 versus run 9L (not shown for the sake of brevity).

However, in the case of the downward VDE situation, the effect of a larger neon injection is different from upward VDs in some respects. This is shown in figure 19 and 20. At first, in spite of the generally much stronger vertical instability in downward VDs, with a larger neon injection, there is a significant RE beam formation (${\sim }8.5$ MA in run 12L) as compared with a negligible beam formed with a lower neon injection (run 12). This is due to the fact that with a higher neon injection, the avalanche time scale is not too slow anymore as compared with the time scale over which the scrape-off causes it to decay. In the case of a lower neon injection, as mentioned earlier, before full avalanche growth occurs, the scrape-off effect dominates and causes RE current decay. Nevertheless, unlike in the upward VDs (runs 1L and 9L), in the case of the downward VDE (run 12L), the RE beam decays via scrape-off soon after reaching the plateau phase. One important consequence of this is the overall faster current decay in run 12L vs run 12, which causes the vertical motion time scale to be faster. In summary, a larger neon injection hastens the vertical plasma motion in a downward VDE, while it slows it down (on average over the whole motion) in the case of upward VDs.

In terms of the effect of a second neon injection and neon flushout, the overall evolution has been observed to be qualitatively similar for the large neon first injection scenario as it was for the lower neon first injection cases described earlier. Such an observation holds for the MD cases (run 1L vs run 2L vs run 5L vs run 4L vs run 6L), the upward-VDE cases (run 9L vs run 10L vs run 11L) and downward-VDE cases (run 12L vs run 13L vs run 14L) as well as with a large neon first injection (figures not shown for brevity). Larger RE beam currents with a larger neon first injection also imply a relatively higher fraction of energy extracted by the REs. Also in terms of energy gain and dissipation by the REs, the effect of a second injection and neon flushout was observed to be qualitatively similar to the cases with a lower neon first injection.

5 Summary and outlook

First axisymmetric simulations of ITER disruptions were presented that include REs, vertical plasma motion (VDEs) as well as massive material injections self-consistently. Effects of a post-TQ RE seed profile, sheath boundary conditions and a more comprehensive treatment of the SOL region have not been considered in the present study. Nevertheless, the work provides important insights into the interplay between RE formation, VDE and RE scrape-off, and timed massive material injections/flushout. Firstly, RE formation not only leads to additional cooling of the plasma but also slows down the decay of current during the disruption and, hence, delays the time to final collapse but at roughly similar final currents (confirms the previous qualitative assessments). While different levels of current flattening can lead to differences in the disruption evolution, the effect in general is found to be significantly weaker due to the fast time scale current profile smoothing-out effect of the highly resistive plasma. Due to this, for example, no noticeable difference was found between flattening within $q=2$ and a significant flattening within the LCFS. Given that high-fidelity 3-D TQ predictions for ITER can be challenging to obtain, the reliability of post-TQ predictions that lack the input of a precise current profile has often been considered questionable. However, the weak effect of the J profile observed in our study indicates that in reality, a precise post-TQ current profile might not be as important.

A neon second injection causes faster RE beam decay, but also a correspondingly faster vertical motion to the wall. Effectively the RE current at final termination does not change significantly due to the second injection. While an earlier injection (during avalanche) enhances the avalanche rate, the subsequent evolution remains similar to the plateau injection. On the other hand, neon flushout not only has a far weaker effect on the overall dynamics but its timing does not make much difference either. Furthermore, even with a neon second injection or flushout, the undissipated RE energy is estimated to be ${\sim }70$ MJ (with a lower neon first injection) and ${\sim }100\unicode{x2013}130$ MJ (with a larger neon first injection), which eventually gets deposited via scrape-off and final collapse. The insensitivity of a neon second injection or flushout in this respect is due to the commensurate increase in the energy channeled to REs along with the energy dissipated via collisions, Bremsstrahlung and synchrotron. This means that a second injection of neon or flushout are not effective in reducing the total RE kinetic energy that gets dumped into the wall (of course 3-D instabilities can change this picture).

In general, the behaviour for the upward VDs (both MD and hot VDEs) have been observed to be qualitatively very similar, except for differences caused by the time delay of TQ in the hot VDE (and the consequent size and current density difference). One important difference in our simulations has been the attainment of $q_{95}<2$ in hot VDEs at a lower neon first injection earlier than any significant RE beam forms. This opens up the possibility for a benign situation of instability-induced stochastic plasma causing partial/full loss of remnant REs (primary seed and partly exponentiated REs). However, with a larger neon first injection, attainment of $q_{95}<2$ occurs only after the full RE beam is formed.

We observe that the downward VDE is strongly unstable after the TQ such that the plasma reaches the final collapse significantly earlier and with a significantly larger current (it must be noted that the word ‘unstable’ is hereby used not in the sense of an exponential growth of a perturbation). Faster gradients in the magnetic field below the X-point is the general reason for the faster time scale of downward VDEs. One of the consequences is that the scrape-off does not allow any significant avalanche growth in the case of a low neon first injection (run 12) such that no RE beam is formed. However, with a larger neon first injection, a high-current (multi-MA) RE beam still still forms before promptly decaying via scrape-off.

Larger neon first injection quantities clearly lead to larger final RE beam currents. The fundamental difference (with respect to a large neon first injection) in the overall evolution of the disruption occurs due to the time scale changes in the vertical motion and RE avalanche. Due to this, for upward VDs, we observe a flat RE plateau phase (as compared with being scrape-off dominated) and an overall slowdown of vertical motion at a larger neon first injection. For downward VDEs, this causes a multi-MA RE beam and also a faster vertical motion. A neon first injection quantity does not however change the qualitative evolution of the disruptions with respect to the effect of a neon second injection and neon flushout.

Finally, it must be pointed out that the relative time scale in ITER of neon flushout as compared with the time for the RE beam plateau to final collapse is an important aspect that is often overlooked when extrapolating the benefits of present day benign terminations to ITER. The typical time for RE beam decay observed in our simulations is ${\sim }40\ {\rm ms}$. So in case neon recombination and outward transport of neon neutrals (flushout) in ITER occurs at a comparable time scale, then the expected benign effects of flushout (as observed in present day experiments) might turn out to be negligible or at most marginal. This calls for an urgent and reliable prediction of the recombination and flushout time scales in ITER. Finally, in the scenario where neon flushout does not lead to benign effects, an edge safety factor drop (and associated stochastic and distributed RE loss) might still hold a possibility for a benign termination. The investigation of such benign terminations in ITER with 3-D MHD simulations using JOREK has been performed and is outside the scope of the present paper (Bandaru et al. 2024b).