Interpretative 3D MHD modelling of deuterium SPI into a JET H-mode plasma

The pre-thermal quench (pre-TQ) dynamics of a pure deuterium ( D2 ) shattered pellet injection (SPI) into a 3MA / 7MJ JET H-mode plasma is studied via 3D non-linear MHD modelling with the JOREK code. The interpretative modelling captures the overall evolution of the measured density and radiated power. The simulations also identify the importance of the drifts of ablation plasmoids towards the tokamak low field side (LFS) and the impurities in the background plasma in fragment penetration, assimilation, radiative cooling and MHD activity in D2 SPI experiments. It is found that plasmoid drifts lead to an about 70% reduction of the central line-integrated density (compared to a simulation without drifts) in the JET D2 SPI discharge considered. Impurities that pre-exist before the SPI as well as those from possible impurity influxes related to the SPI are shown to dominate the radiation in the considered discharge. With inputs from JOREK simulations, modelling with the Lagrangian particle-based pellet code PELOTON reproduces the deviation of the SPI fragments in the direction of the major radius as observed by the fast camera. This confirms the role of rocket effects and plasmoid drifts in the considered discharge and reinforces the validity of the JOREK modelling. The limited core density rise due to plasmoid drifts and the strong radiative cooling and MHD activity with impurities (depending on their species and concentration) could limit the effectiveness of LFS D2 SPI in runaway electron avoidance and are worth considering in the design of the ITER disruption mitigation system.

One of the objectives of the ITER DMS is to avoid the generation of a significant RE population [28].To address this, a two-step SPI scheme has been proposed, with a pure hydrogen (isotope) SPI followed by a (mixed-)neon SPI in a second injection [14,[29][30][31].The idea is that the hydrogen SPI may gradually cool the plasma down to a few hundreds of electron volts via dilution, without immediately triggering a TQ.This allows more time for the hot electrons to equilibrate at an intermediate temperature which could suppress the hot-tail RE generation.In this scheme, the secondary impurity SPI is needed to mitigate the thermal and electromagnetic loads via radiation and avoid large localised heat loads on the plasma facing components.Regarding the hydrogen SPI in the first step, experiments on DIII-D [6] and JET [2,9] have demonstrated that the cooling time, defined as the duration between the arrival of the first SPI fragments at the plasma edge and the peak of the plasma current (I p ) spike (as a result of a global reconnection event), is much longer when using deuterium (D 2 ) SPI than impurity SPI.However, as pointed out in the predictions for ITER [14], this would depend on the level of impurities present in the background plasma, which could radiate strongly, enhance MHD activity and shorten the cooling time.Interpretative 3D MHD modelling of existing SPI experiments is needed to better clarify the underlying physics and validate the numerical model used.
An important process to be considered in the study of D 2 SPI is the possible drift of localised high pressure plasmoids formed by the ablated material of frozen SPI fragments (i.e.ablation plasmoids) towards the tokamak low field side (LFS).This is in line with what has been observed in hydrogen isotope pellet fuelling experiments [32], and references therein].For example, ablation plasmoids have been observed to traverse about 20% of the minor radius in about 10 µs on AUG [33].This outward motion of the ablation plasmoids has been attributed to the E × B drift originating from a vertical polarisation induced inside the plasmoids by the ∇B drifts, where E and B are the electric and magnetic fields, respectively [32,34,35].
3D MHD codes (including JOREK) in principle contain relevant physics for modelling plasmoid drifts self-consistently [36].However, limitations on the spatial resolution make it challenging to quantitatively resolve these drifts in 3D MHD modelling from a computational point of view.This is because the typically enlarged SPI ablation sources used in 3D modelling affect the pressure imbalance between the ablation plasmoids and the background plasma, a driving term for plasmoid drifts [32,34].In the cases where plasmoids do not radiate much (such as the D 2 SPI case considered here), the oversized ablation clouds lead to an underestimated pressure drive, whereas in (mixed-)neon SPI cases they could lead to an underestimation of the radiation, thus an overestimated heating of the plasmoids and an overestimated pressure drive.In addition, an oversized plasmoid means that it would need to bend more magnetic field lines while drifting, thus experiencing an overestimated braking effect from the field line tension.It would also encounter stronger braking originating from 3D external currents that connect the top and bottom of the plasmoid [32], thus an underestimated drift.To circumvent this in simulations with the JOREK code, we propose a 'teleportation' model to account for the particle and energy transfers caused by the drifts, as will be detailed in section 3. 1.In this paper, we will present the first interpretative and quantitative 3D non-linear MHD modelling of disruption mitigation via pure D 2 SPI on JET with JOREK, while detailed JOREK validation against (mixed-)neon SPI experiments is referred to [37,38].We will see that the interpretative modelling does not only confirm the role of background impurities in the radiative cooling and MHD activity of D 2 SPI predicted by previous modelling [14], but also identifies the importance of ablation plasmoid drifts in SPI assimilation and core density rise before the TQ with D 2 SPI.This could also be a possible explanation for the density overestimation in simulations of D 2 SPI experiments on DIII-D with the 0D KPRAD code that does not consider plasmoid drifts [39].
The rest of the paper is structured as follows.We introduce the JET D 2 SPI discharge considered in section 2 and describe the numerical model and simulation setup used in the JOREK modelling in section 3. Simulation results are presented in section 4, including the effects of the ablation plasmoid drifts in section 4.1, effects of impurities in the background plasma in section 4.2 and detailed analyses of the modelled MHD activity and radiation property of the considered JET discharge in section 4.3.Section 4.4 presents a direct simulation of the rocket acceleration of SPI fragments (originating from ablation plasmoid drifts) in the considered discharge, using the Lagrangian particle-based pellet code PELOTON and with inputs from JOREK simulations.Section 5 compares synthetic JOREK bolometry with experimental measurements and discusses possible limitations of the JOREK simulation.Conclusions and outlook are summarised in section 6.Effects of the number of SPI fragments on the pre-TQ dynamics of the considered JET discharge are discussed in appendix.

The JET deuterium SPI discharge considered
The JET discharge considered (#96874) is an H-mode plasma with a plasma current I p ≈ 3 MA, on-axis toroidal magnetic field B 0 ≈ 2.9 T, thermal energy W th ≈ 7 MJ, central electron temperature T e0 ≈ 7 keV and central electron density n e0 ≈ 0.85 × 10 20 m −3 before SPI [2,9].A cylindrical D 2 pellet with a diameter of 12.5 mm, effective length to diameter ratio of 1.54 and total number of deuterium atoms of about 1.4 × 10 23 is shattered and injected into the healthy/non-disruptive plasma.The pellet travels along a flight tube and passes through a microwave cavity (MWC) diagnostic before shattering, which provides information on the pellet integrity and mass contained in different pieces in case of pellet fractures [40,41].In JET #96874, the two spikes of the MWC signal shown in figure 1(a) indicate that the pellet broke into two main pieces before shattering.The mass and thus the number of deuterium atoms contained in each piece can be estimated based on the magnitude of the two spikes and in this case the first piece contains about 25% of the total injected material.
The layout of the SPI system and the main diagnostics available during SPI are illustrated in figure 2, including the vertical bolometry (KB5V), horizontal bolometry (KB5H), KL8 and KLDT cameras (short for the KL8-E8WA and KLDT-E5WE fast visible cameras, respectively), line-of-sight 3 (LoS 3) of the polarimetry (KG4), vertical soft x-ray system (KJ4-T) and a Mirnov coil located near the outer midplane on octant 1 (T001).The toroidal location of the High Resolution Thomson Scattering (HRTS) system is also shown, which provides inputs on the electron temperature (T e ) and density (n e ) before the SPI (figure 6).We note that the KL8 and KLDT cameras are located at different sides of the SPI injection port toroidally, among which the KL8 camera has a direct view of the SPI injection plane.Both cameras have wide-angle views of the JET plasma from just below the horizontal midplane.In JET #96874, the KL8 camera was equipped with a D α filter (656.1 nm), while the KLDT camera was set to receive all lights within the visible spectrum (430 − 730 nm).The frame rates of both cameras were set to 18 kHz in this discharge, with an exposure time of 1 µs and 2 µs for the KL8 and KLDT cameras, respectively.This ensures a good time resolution of the events during SPI.More details about the JET SPI system and key diagnostics involved in disruption studies are described in [2,9,10].
The ablation of D 2 SPI fragments can be inferred from the KL8 camera with a D α filter [42]: excited states of deuterium atoms are populated in the ablatants where ionisation takes place and the Balmer alpha transition of the excited states is recorded.As depicted in figure 1(b) and (c), the rises of the sum intensity of the D α emission at 9.039 s and 9.0425 s (3.5 ms apart) represent the arrival of the two groups of SPI fragments at the plasma scrape off layer (SOL), respectively.We use 9.039 s as the reference time t 0 in the rest of the paper for simplicity.The bulk velocity of each pellet piece can be estimated by the travel time (∆t) and corresponding distance (∆L) from the MWC to the SOL, where ∆L ≈ 4.85 m in JET [10].As marked in figure 1(b), ∆t is about 15.9 ms and 18.2 ms for the two pieces, representing a bulk velocity of about 300 m s −1 and 260 m s −1 , respectively.
One evident phenomenon in the KL8 fast camera video of JET #96874 is a strong outward motion of the first group of D 2 fragments, especially towards the end of their ablation, i.e. when the fragments become small/light.To better quantify this motion, a gradient-based technique has been applied to track the inboard edge of the ablation plume detected by the KL8 camera: a set of bins are defined with respect to the major radius (R) or the normalised poloidal magnetic flux (ψ N ), containing a number of KL8 pixels each; a dataset of the maximum D α emission intensity as a function of space (R or ψ N ) and time is created by taking the maximum emission intensity of all the KL8 pixels contained in each bin and for each KL8 camera frame (i.e. each time slice); the inboard edge of the D α brightness distribution is then located by finding the position that has the maximum gradient of D α emission intensity in the spatial direction.Corresponding results for JET #96874 are shown by the filled circles in figure 3 (the dashed cyan trace is a polynomial fit), where the outward motion of the SPI plume  along R is evident.Note that no correction of possible toroidal expansion or rotation of the SPI plume has been applied in this analysis.This is however not expected to strongly affect the results shown here since no evident toroidal motion of SPI plumes was seen in JET #96874.
The outward motion of the SPI fragments shown in figure 3 is attributed to a 'rocket effect' originating from the drifts of the ablation plasmoids, as will be confirmed in section 4.4 via simulations with a Lagrangian particle-based pellet code.Similar phenomena have been observed in D 2 fuelling pellet experiments on AUG [33].The mechanism of the rocket effect is as follows: drifts of the ablation plasmoids towards the tokamak LFS lead to a stronger shielding and thus reduced heat fluxes reaching the fragment surface on the LFS, causing weaker ablation than on the high field side (HFS) and thus a net force and fragment acceleration along R. The more evident fragment motion towards the end of the ablation in JET #96874 could be explained by the decreasing fragment inertia along with its ablation [43].This could also be a possible explanation for the less evident radial motion of the second group of fragments in this discharge (not shown here for conciseness), which have much larger size thus inertia than the first group (figure 8).The lower plasma temperature and thus weaker ablation experienced by the second group of fragments, due to the initial cooling caused by the first group, could also play a role.
In terms of the pre-existing impurities that could dominate the radiative cooling in D 2 SPI discharges, tungsten (W), neon (Ne), carbon (C), beryllium (Be) and nickel (Ni) were present in JET #96874 according to the charge exchange (CX) and vacuum ultra violet (VUV) spectroscopy measurements (only available before SPI).The existence of Ne could be explained by a series of (mixed-)Ne SPI experiments right before this discharge.Ne and C concentrations are on the order of 1 × 10 17 m −3 based on the CX measurements.Without direct measurements of the other impurity concentrations in this discharge, simulations with the integrated modelling suite JINTRAC [44] have been performed, utilising in particular the coupling between the 1.5D transport code JETTO and the impurity transport code SANCO.With the assumption of homogeneous impurity distributions, scans of the impurity mixture with JINTRAC show that Ne and C contribute mainly to the effective charge (Z eff ) while W dominates the radiation before SPI.Considering only W, Ne and C in the impurity mixture and using a fixed Ne and C concentration of 1 × 10 17 m −3 each, W concentration is estimated to be about 4.5 × 10 15 m −3 to match the measured radiated power (P rad ) and Z eff before SPI.In this case, radiation from pre-existing impurities contributes to about 95% of the P rad before SPI, among which W radiation takes about 80%.
Timing of various events during SPI in JET #96874 is illustrated in figure 4, with a detailed description of each signal given in the figure caption.Line x marks the time when the first group of fragments arrive at the plasma edge, i.e. at t 0 = 9.039 s (figure 1).Lacking a direct measurement of local T e during SPI in this discharge, the line-integrated soft x-ray emission from a LoS that passes through the plasma center (data validated up to about 6.1 ms, blue trace in (a)) is used as an indicator of the core T e evolution.It can be seen that core T e starts to drop from t ≈ 2 ms (line y) along with the penetration and assimilation of the first group of fragments.P rad starts to increase rapidly from t ≈ 6 ms (line z), based on either an inversion using the vertical bolometry and assuming axisymmetric radiation (solid blue trace in (b)) or a 3D tomographic reconstruction without assuming axisymmetric radiation (dashed red) using the Emis3D tool [45,46].W th gradually drops mainly via conduction up to t ≈ 6 ms (line z), after which radiation plays an increasingly important role, as seen from panel (d).A strong burst of MHD activity occurs at t ≈ 6 ms, correlated with the fast increase of P rad , rapid drop of W th and growth of an n = 1 locked mode (n is the toroidal mode number).This marks the onset of the TQ.The I p spike did not occur until t ≈ 10.5 ms, along with another burst of MHD activity and a large locked mode (line {).This marks the end of the TQ in this discharge.In this study we will focus on the pre-TQ dynamics and the TQ onset, i.e. up to the rapid rise of P rad and drop of W th (line z).Note that no RE beams were observed in this discharge, though not shown here for conciseness.and the reconstruction with Emis3D (dashed red), (c) thermal energy evolution derived from diamagnetic loops, (d) thermal energy loss (blue), radiated energy (red), and 20% of the thermal energy loss since t 0 excluding the Ohmic heating contribution (dashed green) to help guiding the eyes, (e) poloidal magnetic field perturbations measured by a fast Mirnov coil and (f ) plasma current and n = 1 locked mode amplitude.

The JOREK model used
The basic JOREK model used comprises ansatz-based reduced MHD equations including parallel (with respect to the direction of the magnetic field) flows and an extension for neutrals and SPI, where the poloidal magnetic flux, toroidal current density, poloidal flow potential, toroidal vorticity, plasma mass density, total plasma (ion and electron) pressure, parallel velocity and neutral mass density are evolved in time, as detailed in [12,16,21].A fixed-boundary setup is used for all harmonics, corresponding to an ideal wall at the computational boundary.The parallel velocity at the intersection of the open field lines with the boundary is fixed to the ion sound speed.The ion temperature (T i ) is assumed to be the same as T e in these simulations.
In terms of the SPI, solid fragments are assumed to travel along their injection trajectories with their initialised velocities.The ablation rate of each fragment (∂N/∂t), defined as the number of ablated deuterium atoms per second is estimated by a neutral gas shielding (NGS) model [47] given as Radiative cooling rates of tungsten (red), neon (blue), carbon (green), deuterium line radiation (dash-dotted magenta) as well as recombination and bremsstrahlung radiation (dashed cyan).
∂N/∂t = 4.12 × 10 16 r 1.33  p n 0.33 e T 1.64 e , where r p is the fragment radius in meters, n e the electron density in m −3 at the fragment's location and T e the local electron temperature in electron volts.The ablated neutrals are then deposited around each fragment with a Gaussian shape in the toroidal, poloidal and radial directions, the half e −1 widths of which are set to 0.5 rad, 15 cm and 2 cm in this study, respectively.This deposition is artificially elongated in the toroidal direction due to limited toroidal resolution that can be reached, as discussed in [14].The sum of all ablated neutrals acts as a source term for the diffusive neutral equation.
Radiation is included as a sink term in the pressure equation.The total radiated power consists of deuterium line radiation (P L ), recombination and bremsstrahlung radiation (P RB ) as well as impurity radiation (P imp ) , i.e.
where n N is the neutral density, n imp,i the density of a given impurity species and N imp the total number of impurity species considered.n e and n N are evolved with time in the JOREK model used, whereas n imp,i is fixed in each simulation and assumed to be homogeneous.L rad,L , L rad,RB and L imp,i represent the cooling rates of the deuterium line radiation, recombination and bremsstrahlung radiation as well as impurity radiation, respectively, with relevant intensities depicted in figure 5. L rad,L and L rad,RB shown by the dash-dotted magenta and dashed cyan traces are directly taken from the open ADAS data [48].We can see that L rad,RB ≪ L rad,L for T e > 2 eV and that L rad,L is comparable to L imp .The contribution to P rad also depends on n e and n N (equation ( 2)) that vary in time during SPI and will be quantified in the simulations (section 4).
In terms of the impurities, L imp of each impurity species (Ne, W and C) shown in figure 5 is estimated using the radiative cooling coefficients of each charge state taken from open ADAS and assuming that the charge state distribution is in coronal equilibrium.W radiates strongly in a wide range of T e , while Ne radiation increases sharply at relevant post-dilution temperature range (around 100 eV [12]).We will thus consider Ne and W as the main impurities in the background plasma in the modelling of JET #96874 in the rest of the paper, while C is shown not to affect the radiation strongly in this discharge, at least at the onset of the TQ (section 4.2).It is worth mentioning that the uncertainties of the impurity transport and distribution as well as the negligence of possible impurity influxes during disruptions [49] allow margins for the effective Ne concentration to be used in JOREK simulations, as will be discussed in section 4.2.
Regarding the drifts of ablation plasmoids, as discussed in section 1, here we propose an ad hoc teleportation model to account for the particle and energy transfers caused by the drifts, given the limitations on the attainable spatial resolution in 3D MHD modelling.This was inspired by a 'backaveraged' model in which the ablated material from pellets is uniformly added to the plasma on the flux surfaces exterior to that where the instantaneous ablation occurs [50].In the teleportation model, the ablated neutrals from each SPI fragment at each time step are shifted from the current location of each fragment (i.e. the birth place of the neutrals) by a certain distance (∆ drift ) along the outward R direction, together with a transfer of energy from each fragment's location to the deposition location of corresponding neutrals.The energy transfer rate corresponds to the total ablation neutral source times a few tens of electron volts (∆E).This is to consider the overall effects of dissociation, ionisation, radiation and early heating before plasmoid drifts, which would have taken place at the fragment's location without teleportation.The setup of ∆ drift and ∆E will be discussed in section 4.1 along with the simulations.
It is worth clarifying that despite the inclusion of the ablation plasmoid drifts via the teleportation model above, the alteration of the fragment trajectory, e.g.via the rocket effect, is not considered in the JOREK simulations shown in this paper for simplicity.This means that the solid fragments are still assumed to travel straight along their injection trajectories with their initial velocities in JOREK.This choice is partially justified by the fact that the deviation from the fragments' ballistic trajectories is not evident for the second group of fragments in the discharge considered, while for the first group, as discussed in section 2, it is not evident either until the fragments are almost fully ablated.The assumption on the fragment motion is thus expected to have limited effects on the JOREK modelling results in this paper.

Simulation setup
The initial equilibrium used in the JOREK modelling is based on a pressure-constrained equilibrium reconstruction of JET #96874 shortly before SPI.The initial flux-surface-averaged n e , T e and q profiles versus ψ N are shown in figure 6 (red  traces), where q is the safety factor.We can see that the input n e and T (assuming T i = T e ) represent reasonably the HRTS and CX measurements and that q < 1 on the magnetic axis.The computational grid on a poloidal plane is illustrated in figure 7, containing 120 radial grid points and 128 poloidal points within the last closed flux surface (LCFS) and about 4000 elements in the SOL.Toroidal harmonics up to n = 10 are included in 3D simulations, while only the n = 0 component is included in axisymmetric/2D simulations (i.e. in appendix).
Plasma resistivity (η) is taken from the Spitzer formulation −3/2 [51] for T e up to 1.7 keV, above which η is set to be uniform.This gives an initial central resistivity η 0 ≈ 5η sp .The parallel heat diffusivity (κ ∥ ) follows the Spitzer-Härm formulation and has a dependence of κ ∥ ∝ T 5/2 e [51], while the perpendicular heat diffusivity is set to 2 m 2 s −1 .The parallel particle transport is purely convective except for a short period around the time when the first post-teleportation plasmoids arrive at the plasma, where parallel particle diffusion is also included for numerical reasons.
The perpendicular particle diffusivity is set to 2 m 2 s −1 .The initial central kinematic viscosity acting on the perpendicular flow is set to 4 m 2 s −1 , while the initial viscosity acting on the parallel flow is about 40 m 2 s −1 .
In terms of the SPI, the bulk velocities of the two groups of fragments are set to 300 and 260 m s −1 , respectively, following the experimental observations discussed in section 2. A ±40% spread is included in the initialisation of the fragment velocities, considering the evolution of the SPI plume length in the experiment as well as the statistics from laboratory measurements [52].The total number of deuterium atoms injected into the plasma is set to 1.4 × 10 23 in the simulations, corresponding to the quantity contained in the original frozen pellet, as specified in section 2. This means that no material loss during the shattering process is considered in this study, which could lead to an overestimation of the injected materials.This choice is due to the difficulty in quantifying the actual injected quantities into the plasma in the experiment, especially when considering the uncertainties of the gas generated during the shattering process [53] and the penetration of gas into the plasma.In fact, the latter is still under active study [54].As will be detailed in section 4, however, plasmoid drifts are expected to play a more important role in the material loss in the D 2 SPI discharge considered here and the negligence of other loss mechanisms is not expected to affect the main results.The number of SPI fragments generated by the shattering process can be estimated by a fragmentation model [53] and is set to 280 and 32 in the two groups, respectively.The effects of the fragment number and sizes have been studied via axisymmetric JOREK simulations (appendix), showing that they do not qualitatively affect the pre-TQ dynamics of the D 2 SPI discharge considered.A histogram summarising the number of fragments at each range of fragment radii is shown in figure 8 (shaded areas), sampled based on the probability distribution from the fragmentation model (dashed red curve).The fragments are initialised such that only small fragments are present in the front and rear of the SPI plume (containing less than 1% of the total mass), in accordance with laboratory measurements [52].

Effects of ablation plasmoid drifts on SPI fragment penetration and assimilation
As discussed in section 1, drifts of the ablation plasmoids could limit the effective assimilation of D 2 SPI fragments injected from the LFS.In this section, we will study relevant effects in the considered JET discharge using the teleportation model introduced in section 3.1.The main parameters in the teleportation model are the teleported energy per atom ∆E and the drift distance ∆ drift , describing the energy and particle transfers caused by the plasmoid drifts, respectively.∆E can be estimated based on previous studies on hydrogen fuelling pellets with dedicated pellet ablation and homogenisation codes such as the hydrogen pellet injection (HPI2) code [43].Specifically, it has been shown that an energy of 20 ∼ 25 eV per atom is typically consumed in the dissociation, ionisation and radiation of D 2 fuelling pellets and that the ablation plasmoids typically reach a few tens of electron volts during drifts (which last on the order of 10 µs) [32,33].This gives an estimate of ∆E about 50 eV.3D JOREK simulations varying ∆E from 50 to 200 eV with a fixed ∆ drift have shown limited effects on the pre-TQ dynamics such as the n el , T e and P rad evolution in JET #96874, though not detailed here for conciseness.We will thus use a fixed ∆E = 50 eV in the rest of this paper.
As for ∆ drift , a scaling law exists for hydrogen isotope fuelling pellets [55], showing a dependence on the pellet velocity, radius, central n e and T e , etc.However, this may not be directly applied to SPI cases considering the modifications of plasma profiles and interactions among different plasmoids that were not included in the scaling law.Given these uncertainties, we estimate the effective ∆ drift in JET #96874 based on the measured central line-integrated electron density (n el ), assuming the same ∆ drift for all ablation plasmoids at this stage.Corresponding simulations with a fixed ∆E = 50 eV, uniform Ne concentration of 1 × 10 17 m −3 and W concentration of 4.5 × 10 15 m −3 are depicted in figure 9.The case with ∆ drift = 0 cm, i.e. no teleportation (red curve in (a)) clearly overestimates the measured n el (blue) in JET #96874, highlighting the important role of plasmoid drifts in this discharge.The simulation with ∆ drift = 40 cm (green trace) captures better the level of measured n el .We note that a simulation with ∆ drift = 0 in the teleportation model here effectively represents a simulation without plasmoid drifts due to the limited spatial resolution discussed before, so we will refer cases with ∆ drift = 0 as cases without plasmoid drifts in the rest of the paper.
The total ablation rate (∂N/∂t) of all fragments is shown in figure 9(b).We can see that ∂N/∂t is much higher for the case with drifts than without drifts before the entry of post-drift ablation plasmoids at t ≈ 0.7 ms.This is attributed to lower dilution and thus higher local T e when considering plasmoid drifts, as can be seen from the n e and T e profiles in figure 10.The solid traces there are from the simulation without drifts, while the dashed cyan and green lines are from the one with ∆ drift = 40 cm.This delayed and much reduced assimilation with plasmoid drifts leads to a higher local T e and thus a much higher ∂N/∂t given its strong T e dependence (equation ( 1)).The ∂N/∂t of the case with ∆ drift = 40 cm drops after the  entry of post-drift ablation plasmoids and remains lower than the case without drifts until t ≈ 4 ms due to smaller fragments remaining after the stronger ablation in the initial phase.
The second group of fragments arrives at the plasma at t ≈ 3.5 ms (section 2) and starts to ablate strongly from t ≈ 4 ms.The difference in the ablation rate is also expected to affect  the effective penetration of SPI fragments.As illustrated in figure 11, the case with plasmoid drifts exhibits shallower fragment penetration than the case without drifts before being fully ablated.
The strong dilution cooling in the case without plasmoid drifts leads to strong MHD activity and an early premature TQ with the first group of fragments, as shown in figure 12 and the red curve in figure 14(b).These are also reflected by the substantial distortion of the electron pressure (P e ) profiles and the collapse of the central T e starting from t = 1.832 ms in figure 10 (solid lines).The premature TQ leads to a loss of about 15% of the original W th , mainly via conduction along stochastic filed lines.Note that P rad remains below 100 MW in this case (red curve in figure 14(a)), the level of which is not expected to dominate the thermal energy loss, as inferred from figure 4(d).JOREK simulations assuming no impurities (thus negligible P rad ) have also been performed to confirm this, though not detailed here for conciseness.
Magnetic energies of the n = 1 − 10 perturbations in the simulation without plasmoid drifts (red case in figure 9 and solid lines in figure 10) are depicted in figure 12. Poincaré plots at a number of selected time slices (marked by the vertical dotted black lines in figure 12) are shown in figure 13.It can be seen that there is a strong increase of all of the harmonics up to t = 0.496 ms, when the vanguard fragments reach the q = 3 surface (figure 13(a)).The mode spectrum becomes broader and a premature TQ is triggered in this simulation at t = 1.832 ms.This in turn leads to lower T e and slower ablation of the second group of fragments than the first group, as seen from the red curve in figure 9(b).The MHD activity is weaker with the second group of fragments (figure 12), with the reformation of magnetic flux surfaces and a dominant m/n = 3/2 mode as seen in figures 13(e) and (f ), where m is the poloidal mode number.This is different from what was observed in the experiment (section 2), where the strong MHD activity and the onset of the TQ did not happen until t ≈ 6 ms (figure 4), i.e. with the second group of fragments.This once again emphasises the importance of plasmoid drifts in this D 2 SPI discharge.We also note here that ∆ drift = 40 cm obtained in the simulations is expected to be a reasonable approximation of the level of plasmoid drifts in #96874.This is based on more recent D 2 SPI experiments on JET, where measurements of the electron density profile during fragment ablation were made possible (not at the time of #96874).By comparing the ablation location (based on the D α filtered KL8 camera) with the material assimilation location (based on the density profile), one can infer the drift distance in a given discharge.Detailed analyses and modelling of the new SPI experiments are planned and will be presented at a later stage.

Effects of impurities on radiative cooling
As discussed in section 2, background W, with a concentration (n W ) of about 4.5 × 10 15 m −3 , dominates P rad before SPI in JET #96874.The contribution from background Ne, with a concentration (n NE ) on the order of 1 × 10 17 m −3 , is expected to play a more important role along the cooling with SPI given its ability to radiate strongly at T e < 150 eV (figure 5).Note that T e < 150 eV can be reached via dilution cooling, especially at the plasma edge, as seen in figure 10(b).Considering the uncertainties of impurity transport and possible impurity influxes from the wall during SPI, thus the uncertainties of impurity concentration and spatial distribution, we have performed scans of n NE (and n C ) with a fixed n W = 4.5 × 10 15 m −3 to determine the effective n NE to match the measured P rad , as summarised in figure 14.All of the simulations except the red case are with ∆ drift = 40 cm, which has proven to match reasonably the measured density evolution (figure 9(a)).Cases with the same ∆ drift but different impurity concentrations exhibit very similar n el , not shown here for conciseness.
We can see from figure 14(b) that Ne plays a key role in the post-SPI radiative cooling in this discharge.In particular, the simulation with n NE = 1 × 10 18 m −3 (cyan trace) exhibits  12).Blue circles: locations of the remaining fragments at each time, with the radii set to five times the fragment radii for illustration purposes.a rapid increase of P rad and fast drop of plasma thermal energy at t ≈ 7 ms, which are absent in the cases with lower Ne concentration.This rapid increase of P rad is correlated with a burst of MHD activities and will be detailed in section 4.3.A possible explanation for the much higher n NE needed to match the measured P rad (versus n NE ≈ 1 × 10 17 m −3 based on pre-SPI measurements) is a rapid impurity influx from the first wall during disruptions [49], especially when considering the loss of the ablation plasmoids (via drifts) in JET #96874 that could cause large heat fluxes onto the wall and a release of impurities via evaporation or sputtering, though it is difficult to quantify this effect with available measurements in this discharge.Note that the impurity influx could contain various impurity species, whereas here we have only varied n NE to represent the overall effect for simplicity.Impurity transport and thus the modification of its spatial distribution, especially in the existence of MHD activities [56], could also add to the uncertainties of the effective post-SPI Ne concentration needed to match the measured P rad .

Modelled MHD activity and radiation property of the considered JET discharge
In this section, we focus on the JOREK simulation with ∆ drift = 40 cm, n NE = 1 × 10 18 m −3 and n W = 4.5 × 10 15 m −3 , i.e. the case that matches best the available density and P rad measurements in JET #96874 (the cyan case in figure 14).The magnetic energy evolution of the n = 1 − 10 perturbations in this simulation is shown in figure 15.We recall that the two groups of solid SPI fragments enter the plasma at t = 0 and t ≈ 3.5 ms in the simulation (section 2), respectively.With a ∆ drift = 40 cm in the teleportation model, the ablated neutrals and thus ionised ablation plasmoids do not enter the computational region (figure 7) until the postablation/remaining fragments penetrate deep enough into the plasma.As indicated by the rises of magnetic energies in figure 15, the post-teleportation plasmoids of the two groups of fragments start to enter the plasma at about 0.7 ms and 4.8 ms, respectively.The MHD activity is stronger with the second group of plasmoids, including a fast growth of n = 1 activity from t ≈ 6 ms, similar to the experiment (figure 4).Higher-n modes grow consecutively, for example n = 2 from 6.5 ms, n = 3 from 6.8 ms and modes with higher n numbers from 7.1 ms.The n = 1 activity remains the strongest in entire phase.
The existence of strong n = 1 activity around the onset of TQ in the simulation is also supported by the helical emission structures recorded by the KL8 and KLDT cameras in JET #96874, as shown in figure 16.The helical blue trace there is obtained by fitting the dominant helical emission structure with the equilibrium magnetic field lines at a given time using the FLARE (Field Line identificAtion fRom imagEs) tool on JET, while the blue dots in the left image are the reference points used for the fitting.The magnetic field line shown by the blue trace is identified to be on the q = 3 surface, indicating a strong cooling of a flux tube with q = 3. Comparing the images from the KLDT camera without any filters (left) and the KL8 camera with a D α filter (right), we can see that D α emission dominates the emissivity within the visible range at this time.We note that other D α emission structures are also visible on the top right of the KL8 image, i.e. near the SPI port.These are expected to originate from SPI fragments at different locations.
Flux-surface-averaged profiles of the JOREK simulation shown in figure 15 are depicted in figures 17 and 19(b)-(f ), focusing on the evolution with the first group and second group of fragments, respectively.It can be seen from figure 17 that n e progressively increases from the plasma edge to the core along with fragment penetration and ablation.The location of the dominantly ablating fragment at each time is marked by a cross symbol with corresponding color.We note here the more outward density peak (thus the assimilation location) with respect to the fragment location (thus the ablation Visible emission structures at 9.04504 s in JET #96874 (i.e.t = 6.04 ms in this paper), taken both from the KLDT camera without any filters (left) and the KL8 camera with a Dα filter (right).The blue trace results from fitting the blue dots shown in the left image with the equilibrium magnetic field lines and is identified to reside on the q = 3 surface.location), as a result of the teleportation/plasmoid drifts.T e gradually decreases mainly due to dilution in this phase, as indicated by the weakly perturbed P e profiles.The radiated power density peaks around the q = 3 surface from t = 2 ms, where T e is cooled down to about 30 eV and the radiative cooling rate of Ne reaches its maximum (figure 5).
As shown in figure 17(e), local current perturbations occur at several main rational surfaces, together with a global current contraction in the outer region (around q = 3).These indicate that both the helical cooling and n = 0 current contraction play a role in the onset of MHD activities here.At t = 2 ms (green line in figure 17(e)), for example, there are local current perturbations at the q = 2 and 5/2 surfaces and a strong current contraction in the outer region.The corresponding Poincaré plot shown in figure 18(a) exhibits a small 2/1 island structure and stochastic regions outside.The 2/1 mode then grows and the stochastic region expands, as shown in figure 18(b).The MHD activities become weaker afterwards (figure 17(a)) and the flux surfaces start to rebuild, for example at t = 3.5ms shown in figure 18(c).The first group of fragments are fully ablated shortly after t = 3.5ms.
Flux-surface-averaged profiles with the second group of fragments are shown in figures 19(b)-(f ).The profiles at t = 4 ms are similar to those at 3.5 ms (cyan traces in figure 17), since the post-drift ablation plasmoids from the second group do not enter the plasma until t ≈ 4.8 ms.Along with the penetration and ablation of the fragments, n e increases, T e drops and P rad rises.Similar to those with the first group of fragments and as shown in figure 19(e), j tor profiles exhibit both local perturbations and global contractions.At t = 5.3 ms, for example, local perturbations exist at the 2/1 surface, while corresponding mode structure can be seen in the Poincaré plot shown in figure 18(d).The outer region, where strong current contraction takes place, is stochastic at this time.The n = 1, 2 and 3 modes grow consecutively (figures 15 and 19(a)), leading to the inward expansion of the stochastic region, as shown in figures 18(e) and (f ).The modes grow rapidly from 7.1 ms and almost the entire plasma region becomes stochastic afterwards.The central T e starts to collapse at t = 7.182 ms (figure 19(c)), while the density within q = 1 remains low.
Similar to that of the first group of fragments, radiation peaks at about 30 eV (figure 19(f the radiative cooling rate of Ne reaches its maximum.This can also be seen in figure 20, where 2D contours of n N , n e , T and radiated power density on the SPI injection plane are shown.Another noticeable feature of the radiation is that the deuterium line radiation (i.e.P L in equation ( 2)) contributes strongly at the later phase of plasma cooling.At t = 7.184 ms, for example, the radiated power density is dominated by deuterium line radiation on the SPI injection plane, as illustrated in figure 21.To have an idea of the integrated quantities considering the entire torus, the time evolution of various radiated power components defined in equation ( 2) is shown in figure 22.We can see that the radiation from impurities (Ne and W in this case) dominates the  15, where time in milliseconds is listed in the legend of (b) and marked by the vertical dotted lines in (a).The location of the fragment with the highest ablation rate at each time is marked by a cross symbol with the same color.Colors of magnetic energies shown in (a) remains the same as figure 15.Note that Te in (c) is plotted in logarithmic scale.
total P rad up to t = 7.1 ms, whereas P L plays a key role in the rapid increase of P rad afterwards and thus the TQ onset.The sharp increase of P L is attributed to the higher neutral density resulting mainly from the recombination at locally cooled regions (initialised by impurity radiation), as can be seen from the n N contours in figure 20(c).High n N blobs are present at the outer region where T e is cooled down to a few electron volts, in addition to the neutral sources from ablating fragments (black circles).

Direct numerical simulation of the rocket acceleration of SPI fragments
In this section, we report results of fully resolved numerical simulations of the rocket acceleration of SPI fragments in JET #96874 using the Lagrangian particle-based pellet code PELOTON and with inputs from the JOREK simulation presented in the previous section.The PELOTON code (the first version was described in [57] as the LP code), is a three-dimensional, scalable, adaptive, massively-parallel particle code for fully-resolved simulations of the ablation of and SPI fragments by hot and runaway plasma electrons based on the Lagrangian particle method [58].It implements the phase transition (ablation) equations on the pellet surface, the low magnetic Reynolds number MHD equations valid in cold and dense ablation clouds, kinetic models for cloud heating by hot thermal background plasma electrons and volumetric heating by REs, equation of state models atomic processes (dissociation and ionisation), radiation models (for Ne and Ne/D mixtures), and models for the transverse motion of ablation clouds across magnetic field lines, in particular the ∇B drift [59].The pellet rocket acceleration is driven by the ∇B drift that creates asymmetry and non-uniform heating of the ablation cloud [60].As a result, the increased pressure on the HFS compared to the LFS leads to pellet (fragment) acceleration where P HFS (P LFS ) is the pressure at the high (low) field side of the pellet, r p is the pellet radius, and m p is the pellet mass.
For the simulation of the pellet rocket effect, the electrostatic shielding model in PELOTON has been improved by implementing a non-uniform model that computes the electrostatic shielding potential along magnetic field lines [61].As a result, the difference of the HFS and LFS pressure in PELOTON comes from three sources.Smaller shielding length of the cloud at HFS compared to LFS results in (a) larger effective background plasma density penetrating the cloud due to smaller magnitude of the electrostatic shielding potential, and (b) weaker attenuation of hot electrons penetrating the cloud at HFS.The third source (c), gradients of the background plasma density and temperature, can also be included in simulations.In some cases (such as pellet moving through a sharp pedestal), they may significantly change the rocket acceleration.PELOTON's simulations of the rocket acceleration agree with experimental observations at ASDEX-Upgrade [62].
For evaluation of the rocket acceleration in JET #96874, we considered a typical SPI fragment injected along the axis of the expected SPI cone (i.e. at 25 • angle to the vertical line).
The background plasma temperature and density experienced by this fragment during the ablation process are taken from the JOREK simulation described in section 4.3 and plotted in figure 23(a).Figure 23(b) depicts the evolution of the fragment radius and the ablation rate.Using PELOTON, we performed simulations of the ablation and rocket acceleration of this reference fragment.We also assumed that an identical fragment,     experiencing the same plasma temperature and density, was injected at 15 • angle to the vertical line as a comparison.
Figure 24 shows the computed trajectories of both fragments due to ablation and rocket acceleration.It can be seen that the computed trajectories and experimental measurements (figure 3) occupy approximately the same spatial domain.More importantly, they exhibit similar deviations in the direction of the major radius.We thus conclude that the rocket acceleration of fragments is the major factor causing deviations of their trajectories in JET #96874.More precise and quantitative comparison of the simulation results with the experimental data is limited by the fact that experimental points represent the front of the entire ablation plume based on the gradient of maximum D α emission intensity observed in the fast camera.Different fragments may dominate the D α emission at different phases of the ablation process, whereas in PELOTON simulations each fragment is tracked separately.

Discussions
Synthetic measurements of the horizontal (KB5H) vertical (KB5V) bolometers for the JOREK simulation presented in section 4.3 are compared in detail with corresponding JET bolometer measurements in this section.The lines of sight of the KB5H and KB5V bolometers are illustrated in figure 25, while the line-integrated radiated power density from the experiment and the JOREK simulation is shown in figure 26.The color scale is kept the same in all the plots in figure 26 for the comparison.We can see that although the simulation captures the overall time evolution, a detailed match with the measured radiation pattern is missing.In particular, as shown in figures 26(a) and (b), the experimental pattern features a more peaked radiation with one maximum at around 160 • on the KB5H plane and two maxima at around 270 • and 245 • on the KB5V plane from t ≈ 6 ms, whereas in the simulation the radiation is more evenly distributed on these planes (also seen in figure 25).
The existence of the peaked radiation in JET #96874 is anticipated to be correlated with the MHD activity during TQ.In [23], similar radiation pattern was observed in a JET argon MGI discharge and found to be related to a large 2/1 activity, qualitatively reproduced by JOREK simulations with the argon source moved in an ad hoc way from the plasma edge into the 2/1 island during the simulation.In JET #96874, a strong cooling of a flux tube with q = 2 is identified at t = 6.4 ms based on the fast cameras (figure 27), following the cooling of q = 3 at t = 6.04 ms that has been discussed in section 4.3 (figure 16).The less peaked radiation in the JOREK simulation shown here indicates a less localised cooling.This could be attributed to the assumption of uniform pre-existing impurities and the absence of an impurity influx in these simulations.The toroidally enlarged neutral sources and thus plasmoids used could also play a role.The unknowns on the impurity distribution and influx during TQ and the limitation on the toroidal resolution hinder a further exploration of these effects at this stage.Visible emission structures at 9.04538 s in JET #96874 (i.e.t ≈ 6.4 ms), taken both from the KLDT camera without any filters (left) and the KL8 camera with a Dα filter (right).The blue trace results from fitting the blue dots shown in the left image with the equilibrium magnetic field lines and is identified to reside on the q = 2 surface.

Conclusions and outlook
The pre-TQ dynamics and TQ onset of a pure D 2 SPI into a JET H-mode plasma (#96874) have been studied in detail based on experimental analyses and interpretative 3D MHD modelling with the JOREK code.The interpretative modelling has captured the overall evolution of the post-SPI density and radiated power, although a more detailed match is hindered by the uncertainties of pre-existing and influx impurity distribution and the limitation on toroidal resolution.More importantly, the simulations have highlighted the key role of the ∇Binduced ablation plasmoid drifts towards the tokamak LFS and the existence of impurities in the background plasma in fragment penetration, assimilation, radiative cooling and MHD activity in D 2 SPI experiments.
The limitation on spatial resolution has motivated us to develop and implement an ad hoc teleporation model in JOREK to consider the particle and energy transfers caused by the drifts of ablation plasmoids.By performing and comparing JOREK simulations with and without drifts, it has been demonstrated that plasmoid drifts led to an about 70% reduction of the central line-integrated density in JET #96874.These drifts have been shown to affect the MHD activity and the onset of the TQ in JET #96874 as well.Specifically, plasmoid drifts give rise to lower dilution than the case without drifts during the assimilation of the first group of fragments, leading to higher T e , faster ablation and stronger MHD activity with the second group of fragments and eventually the TQ onset observed in the experiment.Impurities in the background plasma have been shown to dominate the measured radiated power in JET #96874, with tungsten being the main contributor before SPI and neon the key contributor after SPI given its ability to radiate more strongly at lower temperature.Neon concentration has proven to be crucial to recover the measured radiated power evolution and the TQ onset in the considered D 2 SPI discharge.In fact, the strong plasma cooling initialised by neon led to strong recombination, higher neutral density (in addition to fragment ablation), a fast rise of deuterium line radiation and thus a sharp increase of the total radiated power.
Outputs from the JOREK simulation, in particular the local electron temperature and density seen by a given fragment along its penetration have been used as inputs for the Lagrangian particle-based pellet code PELOTON to study the rocket acceleration of SPI fragments in JET #96874.The computed fragment trajectories exhibited similar deviation in the direction of the major radius as observed by the fast camera.This has on one hand confirmed the role of the rocket effects and plasmoid drifts in JET #96874, and on the other hand reinforced the validity of the JOREK simulation results obtained.
The limited core density rise due to plasmoid drifts and the strong radiative cooling and MHD activity in the existence of background impurities (depending on the impurity species and concentration) could limit the effectiveness of LFS D 2 SPI in RE avoidance and are worth considering in the design of the ITER DMS.To better quantify the key parameters at play in plasmoid drifts, a dedicated study using simpler JOREK setup is planned.This will also help understand better the presence of plasmoid drifts in recent JOREK simulations of low Ne/D mixture ratio SPI into ITER H-mode plasmas [63], whereas no evident drifts are observed in the JOREK simulations shown in this paper (without the teleportation model).Factors such as plasma thermal energy and injection geometry (equatorial SPI in the aforementioned ITER simulations versus near top injection on JET) are expected to play a role and will be explored in the dedicated study.
Finally, it is worth mentioning that despite the unfavorable features of the LFS pure D 2 SPI in terms of ablation plasmoid drifts, recent studies with dedicated pellet ablation and homogenisation models have predicted that even a mixture of neon as low as 2% in the pellet could substantially reduce the level of plasmoid drifts, thanks to its much stronger cooling, lower plasmoid pressure and thus a weaker drive for plasmoid drifts [34].This could offer an alternative to the first-step pure hydrogen SPI in the two-step injection scheme proposed for RE avoidance [14,29].Dedicated SPI experiments have been performed on JET recently and 3D MHD simulations are envisaged to fully examine this.

Appendix. Effects of fragment number on pre-thermal quench dynamics
As discussed in section 3.2, the number of SPI fragments generated by the shattering process and the corresponding probability distribution of fragment sizes can be estimated based on a fragmentation model [53], i.e.where d represents the fragment diameter in centimeters, f (d) the relative probability of fragment sizes and K 0 a modified Bessel function of the second kind.α ≡ X R /D, β ≡ X R /(LC) and X R ≡ v 2 perp /v 2 thres , where v perp is the component of the pellet velocity perpendicular to the shattering plate, v thres the perpendicular impact threshold velocity, D the pellet diameter, L the pellet length and C the pellet material constant (2.5 for deuterium).It can be seen that fragment sizes and thus the total number of fragments are very sensitive to the pellet velocity and length for a given pellet species.Considering this sensitivity and the uncertainties in the way the pellet broke (into two main pieces) in JET #96874, i.e. the effective pellet length to be used, we have performed studies on the effects of the number of fragments and corresponding fragment sizes on the pre-TQ dynamics via axisymmetric JOREK simulations.
v thres used in equation (A.1) was found to be smaller than 22.8 m s −1 for deuterium pellets based on laboratory measurements, with a suggested value of 20 m s −1 [53].In terms of the effective length of each piece to be used in the fragmentation model, here we consider two extreme situations of pellet fracture to cover a wide range: type A where the pellet breaks along its length and generates two pieces with the same D = 12.5 mm and type B where the pellet breaks into two pieces with the same length as the original pellet (12.5 mm × 1.54 = 19.25 mm) but different effective diameters.The effective length (for type A) and diameter (type B) of each piece are then estimated based on the mass of each piece, which can be quantified via the total injected mass and the mass ratio given by the microwave cavity diagnostic (section 2).The resulting number of fragments and corresponding characteristic fragment size of each group for the cases involved in the sensitivity study are summarised in table A1.The much larger number of fragments in the first group than the second group, especially when assuming the pellet breaks along its length (cases #1 and #2 in table A1), are attributed to its higher velocity (300 versus 260 m s −1 ) and shorter length used in the fragmentation model.
Results from the 2D scans are depicted in figure A1, including the cases listed in table A1 as well as the 280/32 setup used in the paper.All the simulations have the same setup except the fragment number and sizes, with n NE = 1 × 10 18 m −3 and n W = 4.5 × 10 15 m −3 .Here we only focus on the effects on fragment ablation and assimilation (thus the density evolution) since the highest uncertainty of the fragment number and sizes lies in the first group of fragments (table A1), during which radiation does not dominate (section 4).The choice of 2D modelling can be justified by the fact that 2D and 3D simulations exhibit very similar ablation and n el evolution and that it reduces computational cost substantially.
It can be seen from figure A1(a) that a smaller number of fragments in the first group leads to a higher n el with a fixed ∆ drift = 40 cm.This is because a small number of big fragments has a smaller total surface area thus slower ablation than the cases with a large number of tiny fragments (with the same injection quantities), giving rise to deeper fragment penetration and stronger density rise.A larger ∆ drift would recover better the measured n el in this case.More importantly, as displayed in figure A1(b), simulations without plasmoid drifts strongly overestimate the measured n el regardless of the number of fragments, confirming again the importance of plasmoid drifts in the considered D 2 SPI discharge.We thus conclude from this section that the uncertainties in the number of fragments and fragment sizes only have limited quantitative effects on the modelling of JET #96874 and do not affect the main conclusions of this paper.

Figure 1 .
Figure 1.(a) Microwave cavity signal of JET #96874, (b) sum intensity of Dα emission from the KL8 camera and (c) corresponding KL8 camera images at 9.039 s and 9.0425 s.

Figure 3 .
Figure 3. Tracked inboard edge of the SPI plume based on the maximum gradient of the Dα emission intensity detected by the KL8 camera in JET #96874.Time in milliseconds shown in the colorbar is shifted by t 0 = 9.039 s, i.e. the time when the first group of fragments arrive at the plasma.Outlines of the expected SPI plume from the injection geometry are shown by the red dashed lines.

Figure 4 .
Figure 4. Overview of JET #96874 during SPI, with the time axis shifted by t 0 = 9.039 s.(a) Soft x-ray signal from the LoS passing through the plasma center (blue) and the measured loop voltage (red), (b) radiated power from the vertical bolometry (KB5V, blue)and the reconstruction with Emis3D (dashed red), (c) thermal energy evolution derived from diamagnetic loops, (d) thermal energy loss (blue), radiated energy (red), and 20% of the thermal energy loss since t 0 excluding the Ohmic heating contribution (dashed green) to help guiding the eyes, (e) poloidal magnetic field perturbations measured by a fast Mirnov coil and (f ) plasma current and n = 1 locked mode amplitude.

Figure 6 .
Figure 6.Initial flux-surface-averaged profiles used in the JOREK modelling of JET #96874.Left: ne used in JOREK (red) and from TS measurements at 8.925 s (blue crosses); middle: T = T i = Te used in JOREK (red), Te at 8.925 s from TS (blue crosses) and T i at 8.868 s from CX (green squares); right: q profile used in JOREK.

Figure 7 .
Figure 7. Poloidal computational grid used in the JOREK modelling.

Figure 8 .
Figure 8.A histogram summarising the number of fragments at each range of fragment radii used in the JOREK modelling of JET #96874 (shaded cyan areas) and the probability distribution from the fragmentation model (dashed red curve) for (a) the first group and (b) second group of fragments.

Figure 9 .
Figure 9. (a) Synthetic central line-integrated density, (b) total ablation rate of all fragments and (c) total number of ablated deuterium atoms of simulations without drifts (red) and with ∆ drift = 40 cm (green), respectively.

Figure 10 .
Figure10.Flux-surface-averaged profiles with the first group of fragments, where time in milliseconds is shown in the legend and the location of the fragment with the highest ablation rate at each time is marked by a cross symbol with corresponding color.Solid: simulation without drifts (red case in figure9); dashed green and cyan in (a)-(c): ∆ drift = 40 cm (green case in figure9) at t = 0.978 ms and 3.015 ms, respectively.Note that Te in (b) is plotted in logarithmic scale.

Figure 11 .
Figure 11.Averaged ψ N location of the remaining fragments from the first group for the simulation without (red) and with plasmoid drifts (green), corresponding to the cases shown in figure 9.

Figure 12 .
Figure 12.Magnetic energies of the n = 1 − 10 perturbations for the simulation without plasmoid drifts (red case in figure9).The time slices marked by the vertical black lines correspond to those shown in figures 10 (up to t = 1.832 ms) and 13.

Figure 13 .
Figure13.Poincaré plots at various time slices of the simulation without plasmoid drifts (figure12).Blue circles: locations of the remaining fragments at each time, with the radii set to five times the fragment radii for illustration purposes.

Figure 14 .
Figure 14.(a) Total radiated power and (b) plasma thermal energy evolution taken from measurements (solid blue traces) and from JOREK modelling scanning Ne and C concentrations with a fixed W concentration of 4.5 × 10 15 m −3 .The red and (dashed) green cases are from the same simulations shown by the red and green traces in figure 9, respectively.

Figure 16 .
Figure 16.Visible emission structures at 9.04504 s in JET #96874 (i.e.t = 6.04 ms in this paper), taken both from the KLDT camera without any filters (left) and the KL8 camera with a Dα filter (right).The blue trace results from fitting the blue dots shown in the left image with the equilibrium magnetic field lines and is identified to reside on the q = 3 surface.

Figure 17 .
Figure 17.Flux-surface-averaged ne, Te, Pe and jtor profiles with the first group of fragments of the simulation shown in figure 15, where time in milliseconds is listed in the legend of (b) and marked by the vertical dotted lines in (a).The location of the fragment with the highest ablation rate at each time is marked by a cross symbol with the same color.Colors of magnetic energies shown in (a) remains the same as figure 15.Note that Te in (c) is plotted in logarithmic scale.

Figure
Figure plots at various time slices of the simulation shown in figure 15.

Figure 19 .
Figure 19.Flux-surface-averaged ne, Te, Pe and jtor profiles with the second group of fragments of the simulation shown in figure 15, with time in milliseconds listed in (b) and marked by the vertical dotted lines in (a).The blue and red traces are taken from the blue and cyan traces in figure 17, respectively.The location of the fragment with the highest ablation rate at each time is marked by a cross symbol with the same color.Colors of the magnetic energies shown in (a) remains the same as used in figure 15.Note that Te in (c) is plotted in logarithmic scale.Profiles at t = 0 are shown here as a reference.

Figure 20 .
Figure 20.From left to right: n N , ne, Te and radiated power density contours on the SPI injection plane at (a) t = 5.3 ms, (b) t = 7.1 ms and (c) t = 7.184 ms.Cyan contours: location of Te = 30 eV; black circles: remaining SPI fragments at each time.

Figure 21 .
Figure 21.(a) Neon, (b) tungsten, (c) deuterium line and (d) total radiated power density on the SPI injection plane at t = 7.184 ms, i.e. the same time as figure 20(c).Cyan contours: location of Te = 30 eV; black circles: remaining SPI fragments at this time.

Figure 22 .
Figure 22.Time evolution of various radiated power components of the JOREK simulation shown in figure 15.The dashed cyan curve is the same as the cyan trace in figure 14(a).

Figure 23 .
Figure 23.(a) Electron temperature and density of the background plasma at instantaneous locations of the fragment as it moves along the injection trajectory, taken from the JOREK simulation discussed in section 4.3.(b) Time evolution of the fragment radius and ablation rate from PELOTON.

Figure 24 .
Figure 24.Computed trajectories of fragments due to ablation and rocket acceleration using PELOTON.Solid line: fragment injected along the axis of the SPI cone, at 25 • angle to the vertical line.Dashed line: the same fragment injected at 15 • to the vertical line.

Figure 25 .
Figure 25.Poloidal cross-section of the radiated power density on the plane of (a) the horizontal (KB5H) and (b) vertical (KB5V) bolometers at t = 6.4 ms of the simulation discussed in section 4.3.The white lines are the lines of sight of JET KB5H and KB5V, with extreme poloidal angles indicated.

Figure 26 .
Figure 26.Experimental (top) and synthetic (bottom) bolometry measurements for the horizontal (KB5H) and vertical (KB5V) bolometers.The color scales represent log 10 the line-integrated radiated power density in (Wm−2).The two vertical dotted lines in (a) and (b) mark t = 6.04 ms and 6.4 ms, i.e. the time slices shown in figures 16 and 27, respectively.

Figure 27 .
Figure 27.Visible emission structures at 9.04538 s in JET #96874 (i.e.t ≈ 6.4 ms), taken both from the KLDT camera without any filters (left) and the KL8 camera with a Dα filter (right).The blue trace results from fitting the blue dots shown in the left image with the equilibrium magnetic field lines and is identified to reside on the q = 2 surface.

Figure A1 .
Figure A1.Central line-integrated density for JOREK 2D simulations (a) with ∆ drift = 40 cm and (b) ∆ drift = 0 cm, i.e. without plasmoid drifts.The number of fragments in the first/second group is shown in the legend.

Table A1 .
List of the number and characteristic sizes of fragments used in the sensitivity study.