MHD modeling of shattered pellet injection in JET

Nonlinear 3D MHD simulations of shattered-pellet injection (SPI) in JET show prototypical SPI-driven disruptions using the M3D-C1 and NIMROD extended-MHD codes. Initially, radiation-driven thermal quenches are accelerated by MHD activity as the pellet crosses rational surfaces, leading to a radiation spike, global stochasticization of the magnetic field, and a complete thermal quench. Eventually, current quenches, preceded by a current spike are seen as the Ohmic heating becomes equal to the radiative cooling. The results are qualitatively similar for both a single monolithic pellet, pencil-beam model, and a realistic shatter to represent the SPI plume. A scan in viscosity from 500 to 2000 m2 s−1 for MHD simulations finds that reducing viscosity increases MHD activity and decreases thermal quench time slightly. A realistic cloud of fragments modeling shows that mixed-D–Ne pellet travels deeper into the plasma core before the thermal quench. At the slow pellet speeds, the pellet is found to be moving slowly enough inward that even the 5% neon in the mixed pellet is enough to effectively radiate the thermal energy available. Radiation toroidal peaking is predicted to be at levels consistent with experimental observations and reduced as the pellet travels deeper into the plasma. These simulations lay the ground work for more-sophisticated validative and predictive modeling of SPI in JET using both M3D-C1 and NIMROD.


Introduction
Due to the large stored magnetic and thermal energy, ITER [1] will need to mitigate the effects of a disruption, which cause a rapid loss of the stored energy due to a plasma instability. The primary candidate for an ITER disruption mitigation system (DMS) is shattered pellet injection (SPI). SPI launches a pellet of high atomic number Z (e.g. neon, argon, etc) material, often mixed with a hydrogen species, to radiate the plasma energy, distributing the energy onto the entire first wall, rather than just the divertor. Before entering the plasma, the pellet is shattered to protect the inner wall which could be damaged by an intact pellet, and simultaneously improves assimilation of particles by generating a high surface area spray of pellet fragments.
Experiments on the largest tokamak JET have been performed to assess the effectiveness of SPI disruption mitigation. These experiments scanned the amount of neon content in an SPI disruption with a variety of deuterium-neon mixtures [2], finding that cooling times and current quench duration are reduced as the impurity content is increased. While synthetic diagnostics were needed to estimate the total levels of radiation due to limited bolometry coverage, they concluded that toroidal peaking factor (TPF = the radiation at the injection location divided by the toroidal mean) had an average value 1.6-1.9 before the current quench and 2.0-2.6 at peak radiation.
Sophisticated modeling is required to assess the efficacy of the ITER DMS, and validation against present-day experiments is necessary to provide confidence in the predictive capabilities of the models. The physics of DMS with SPI involves ablation and radiation physics in addition to conventional plasma dynamics. To model these effects, this work uses the MHD codes M3D-C1 [3,4] and NIMROD [5,6] that have been coupled to a pellet ablation and ionization and radiation subroutines from KPRAD [7]. Benchmarking these codes is critical in improving the models. Modeling with both M3D-C1 and NIMROD allows for verification that the models have been coupled consistently together. The first benchmarking exercise involved a 2D axisymmetric simulation which injected a core impurity source at a constant rate and showed nearly identical results between the two codes [8]. This coupling of MHD to the SPI model has been validated showing good agreement with DIII-D experiments [9]. However, validation with experiments varying major radius and stored energy is needed to ensure accurate SPI predictions for ITER motivating testing of the modeling on JET. While these JET experiments have a large plasma volume of 75 m 3 and stored thermal energy 3-4 MJ, it is still a large extrapolation towards the ITER baseline scenario which has a projected plasma volume of 830 m 3 and stored thermal energy 350 MJ. A similar coupling has also been performed between an SPI model and a reduced MHD model with the JOREK code [10], which found fairly quiet disruptions on JET until the pellet reached the q = 1 surface, where q is the plasma safety factor. Later work using JOREK showed that the simulation of a JET L-mode exhibited a strong toroidal asymmetry in the radiation power density peak [11].
The structure of the paper is as follows. Section 2 shows the setup of MHD simulations for SPI on JET, and compares the differences between NIMROD and M3D-C1 implementations. Initial SPI disruption simulations using a monolithic and pencil beam representation find typical SPI disruptions to have a modest sensitivity to viscosity, and are discussed in section 3. In section 4, SPI simulations with a realistic plume representation are performed for pure-Ne pellet and mixed-D-Ne pellet, and with the finding that the mixed-D-Ne pellet travels deeper into the core before the end of the thermal quench. Additionally, the thermal quench times become similar between the two pellets at lower velocities, and radiation toroidal peaking is reduced as the pellet travels deeper into the plasma. Conclusions are given in section 5.

Simulation setup
In this work, we use impurity-modified single-fluid nonlinear 3D resistive MHD equations coupled to ionization and radiation subroutines from KPRAD [7] and a pellet ablation model [12]. The density continuity equation is modified to include the source (S) from the ablated pellet. Then the MHD momentum equation and Ohm's law are the standard MHD equations but are modified indirectly by the ablated pellet which influences the resistivity, and the mass density which is summed over the charged electrons, charge ions, and neutrals. The NIMROD [6,9] and M3D-C1 [4] implementations differ in the discretization and temperature equations. NIMROD is discretized with finite elements in the poloidal plane and Fourier decomposed in the toroidal direction, whereas M3D-C1 uses finite elements with C1 continuity in all dimensions. To calculate the temperatures, NIMROD advances a single energy equation summed over all charged species which assumes instant thermalization: where the temperature T is the same for all species (T e = T D = T Ne = T). The subscripts e, D, and Ne correspond to electrons, deuterium, and neon, respectively. Here n tot = n e + n D + n Ne is density summed over all species including neutrals. The variable P rad is the radiation power density calculated by KPRAD, V is the single fluid velocity for charged, neutral, and electrons, Γ = 5/3 is the adiabatic index, q is the heat flux (which includes anisotropic parallel and perpendicular thermal conduction). The heat source Q ion includes ionization and recombination which are sources and sinks for charge states calculated by KPRAD.
M3D-C1 uses a two-temperature model with the electron energy equation: and the ion energy equation: Here n i = n D + n Ne is the summed ion density, Q ei is the electron-ion exchange energy, Π is the stress tensor, and the 1 2ω V 2 term accounts for the net loss of kinetic energy caused by the slowing of the fluid velocity as new particles are added.
In absence of strong parallel flows which are not observed in this work, the viscous heat is small. So the largest difference between the two models is two temperatures versus one. With the two-temperature model, the radiated energy from the neon is subtracted from the electron temperature equation, with the ions cooled by collisional exchange. This could potentially slow the time of the thermal quench compared to the instantaneous thermalization of NIMROD. However, T e gets very cold in the presence of a strong radiative source, which significantly increases the electron-ion collision frequency and Q ei . Thus, the difference in models is unlikely to have a large effect. Note that while two MHD codes are used, this paper is not meant to be a detailed benchmark that requires careful setup, and thus simulation parameters between the two codes will vary. While this makes it more difficult to assess differences, similarities between the codes can be considered more robust as the result does not depend on the specific simulation setup.
We treat the shatter pellet fragmentation with multiple levels of fidelity. At the simplest level, the shatter is ignored and the pellet is treated as one single monolithic sphere. As the next level of complexity, the pellet can be treated as a pencil beam of pellet fragments where the fragments enter the plasma in a straight line trajectory with a fixed delay between fragments. For the highest level of complexity, a fracture threshold shattering model is used to create realistic SPI plumes based on the size, speed, and composition of the pellet, as well as the geometry of the SPI system [13]. In this paper, a simplified version of this model is used with a uniform pellet radius, calculated such that the full shattered cloud have the same total ablation rate at fixed density and temperature. This simplification allows for a reduced number of fragments to represent the plume. Trajectories of the shattered pellet launched from JET barrel B are shown in figure 1, where each pellet fragment is deposited as a neutral density in the form of a Gaussian on the poloidal plane and a Von Mises distribution toroidally.. All pellet trajectories in this paper are injected from JET barrel B. Kinetic EFIT equilibrium reconstructions have been generated based on target Scenario 1 high thermal energy JET discharge: 95707, which has a plasma current of I p = 2.4 MA. Simulations in this paper will use two equilibria from this discharge at t = 50.875 s with E th = 3.2 MJ and t = 50.544 s with E th = 3.4 MJ. The pressure (p) and safety factor (q) of the equilibria are shown in figure 2 with the q = 2 surface being at ρ = 0.75 and q = 1.5 at ρ ≈ 0.62 for both equilibria. The electron density has an H-mode density profile with pedestal of 0.55 × 10 20 m −3 and an on-axis density of 0.72 × 10 20 m −3 . The poloidal flux contours of the t = 50.544 s equilibrium are shown in figure 1. Note that, while not shown, the contours of t = 50.875 s equilibrium are nearly identical to the t = 50.544 s equilibrium.

Analysis with single and pencil beam
In this section, initial SPI simulations are performed using the single monolithic pellet, pencil beam, and more realistic plume representations, showing that simulation parameters such as viscosity and the plume distribution do not qualitatively change the simulation results. The first 3D nonlinear MHD simulation shows typical characteristics of an SPI disruption. This simulation injects a single monolithic pure-Ne fragment with velocity v f = 150 m s −1 and radius r f = 4.05 mm, viscosity of ν = 1000 m 2 s −1 , and is performed using the M3D-C1 code. Figure 3 shows time-traces of the thermal stored energy (E th ), plasma current (I p ), injected number of electrons (∆N e ), radiated plus ionization power (P rad + P ion ), and Ohmic power (P Ohm ) for the simulation. Additionally, figure 3 shows a time trace of the magnetic energy harmonics for toroidal mode numbers n = 1-8. As the pellet enters into the plasma, the neon atoms ablate, become ionized, and begin to radiate energy. At t = 2.1 ms, the pellet crosses the q = 2 rational surface. There is no significant MHD event until just before t = 4 ms, shown by the peak in perturbed n = 1 magnetic energy. This MHD event precedes a spike in P rad , and a rapid loss in E th . The plume traveling beyond the key rational surfaces may be due to the delay between when the plume reaches the surface and the temperature collapse, as observed in [14]. As the pellet travels farther into the plasma with ρ pellet = 0.42, the remaining thermal energy is terminated by a secondary MHD event at t = 5.5 ms. After the thermal quench, the Ohmic heating becomes equal to the radiative cooling and a predicted slight spike in I p is observed, which is consistent with I p spikes that are commonly observed in SPI disruption experiments on JET.
Like M3D-C1, NIMROD 3D nonlinear MHD simulations injecting a pure-Ne pellet also show a typical SPI disruption. The simulation setup was similar to the M3D-C1 simulation with a faster velocity and slightly larger fragment  where G i is the pellet fragment ablation rate, and i and n are the individual and total pellet fragments. This simulation is qualitatively similar to the M3D-C1 simulation which has a slow initial reduction of thermal energy, followed by two distinct MHD events causing a complete thermal collapse. This simulation has a faster thermal quench than the single fragment M3D-C1 simulation with a thermal quench time of 3 ms, due to the faster fragment velocity. The incoming pellet fragment finds a temperature balance between the fragment traveling deeper into the hot core and the fragment locally cooling the plasma from the radiation, as seen in the middle panel of figure 4. Additionally, the acceleration of the thermal quench can be clearly observed by examining the temperature that the pellet experiences. At t = 2 ms, the local T e of the pellet rapidly drops, suggesting that the plasma is moving away from the pellet. Also shown in figure 4 are NIMROD simulations of SPI injection using the pencil beam and realistic plume representations, finding that increasing the number of fragments and spreading out the fragments causes the thermal quench to be completed at a larger ρ abl . For the pencil beam, the shatter is broken into 50 fragments spread over 10 cm. As multiple fragments are injected into the plasma, each fragment experiences its own ablation from the background plasma temperature, as shown in the middle panel of figure 4. Each fragment balances at a colder temperature than the single fragment simulation with more atoms entering into the simulation, and thermal quench occurring at a larger ρ abl = 0.63. However, some of this difference could be due to the farther inward starting position of the monolithic pellet. All three complexity levels qualitatively exhibit the same effect of the ablation halting before the end of the thermal quench due to the MHD instability and the plasma kinking away from the pellet. Even though the thermal quench occurs at a larger ρ abl for the pencil beam simulation, the thermal quench ends at t = 2.8 ms for both simulations. Like the pencil beam simulation, the realistic plume travels a shorter distance into the core before the completion of the thermal quench. As the realistic plume pellet was injected at a slower velocity which would allow MHD instabilities more time to grow and disrupt the plasma, it is difficult to determine whether the earlier termination is due to the change in the plume representation or the slower pellet velocity.
The previous simulations used a viscosity ν = 1000-2000 m 2 s −1 which was chosen to suppress numerical instabilities. Reducing ν is found to moderately change the thermal quench process by increasing the MHD activity, which increases peak radiation power and decreases thermal quench time. NIMROD SPI simulations with the pencil beam representation are performed scanning three viscosities to assess the importance of this numerical parameter, and is shown in figure 5. Early in time (0-1 ms), P rad and the reduction of E th are similar between all the simulations. As noted in figure 3, the peaks in radiated power correlate with the peaks in MHD amplitude which deposit heat flux into the impurity radiation zone. As the pellet travels deeper into the core and P rad increases at t = 1.5 ms, the simulations begin to differ with the ν = 500 m 2 s −1 simulation having the highest level of radiation and MHD activity. By t = 1.8 ms, the radiation power from the ν = 1000 m 2 s −1 and ν = 2000 m 2 s −1 simulations catches up to the ν = 1000 m 2 s −1 simulation with P rad = 2 GW. At t = 2.1 ms, the radiation levels begin to differ again as the 500 m 2 s −1 simulation's radiation begins to peak. Note that these SPI simulations are numerically expensive towards the end of the thermal quench, especially at lower viscosities. Thus the lower-viscosity simulations in figure 5 were not carried out to the end of the thermal quench. Extrapolating out to the end of the thermal quench, the 500 m 2 s −1 simulation would have a complete thermal quench around t TQ ≈ 2.6 ms. The larger viscosities have slower thermal quenches with an extrapolated t TQ ≈ 2.8 ms for ν = 1000 m 2 s −1 , and an observed t TQ = 2.9 ms for ν = 2000 m 2 s −1 . While qualitatively similar, the thermal quench time is reduced from approximately t = 2.9-2.6 ms as the viscosity is reduced from ν = 2000-500 m 2 s −1 . The lower-viscosity simulations also have stronger radiation peaks.

Simulation with realistic SPI plume
Simulations using both M3D-C1 and NIMROD models indicate that a mixed-D-Ne pellet penetrates deeper into the core than a pure-Ne pellet before the thermal quench is completed. This is due to both the material composition and the higher speed of the mixed pellet. Two types of pellets were simulated using a realistic shattered fragment representation, as illustrated in figure 1. Details of the two pellet parameters have been added to table 1. The first type is a pure-Ne pellet broken into 30 identical fragments with a radius of 1.71 mm and velocity of 150 m s −1 . The second type is a mixed-D-Ne pellet consisting of 5% neon and 95% deuterium by molar mass. To be consistent with JET experiments, the mixed-D-Ne pellet is injected v f = 300 m s −1 . Due to the composition and the faster speed, the pellet is shattered into 85 identical fragments with a radius of r f = 1.21 mm. Both of these pellets have enough neon to rapidly thermal quench the plasma if the neon gets to the core with the pure-Ne pellet having 99.6 × 10 20 neon atoms and the mixed-D-Ne pellet having 18.5 × 10 20 neon atoms. This is compared to the total number of electrons in the initial equilibrium: 46.4 × 10 20 . All fragments in these plumes have the same size and velocity in order to reduce the needed number of simulated fragments, which was determined by analyzing the shattering of the pellet with a fracture-threshold model from Gebhart et al [13]. This gives clouds of fragments with a distribution of sizes. A weighted average was then taken such that the cloud with distributed sizes had the same total ablation rate as the given clouds with uniform fragment size, subject to constant background plasma density and temperature. Based on figure 4, the simplification of the plume is not likely to significantly change the predicted thermal quench. Note that a substantial fraction could be converted to a prefrontal gas which this work does not account for. This prefrontal gas could potentially cool the plasma ahead of the shattered pellet triggering MHD modes and accelerating the thermal quench. To estimate the amount of gas that would be made in the shattering process, the ratio of normal impact kinetic energy to threshold kinetic energy (X r ) is shown in table 1. Estimating figure 5 of Gebhart et al the fraction of the pellet that is turned into gas is between 5% and 20%. Similar JET discharges have been identified with a pure-Ne SPI injection discharge 94579 and a mixed-D-Ne SPI discharge 94575. Figure 6 is the reconstructed P rad and the reconstructed (TPF = P rad (ϕ = 0)/P rad (n = 0)). The TPF is reconstructed using Emis3D [15], assuming that radiation is peaked at the injector location, and falls off as an asymmetric Gaussian. Note that these Emis3D results have significant uncertainties, as they draw from a limited radiation structure pool of helical flux tube shapes. To help quantify the uncertainty, error bars of fit have been included in figure 6, which are estimated by considering the distribution of predicted radiation structures within the fitting pool. The fitted peak radiation is P rad ≈ 1 GW at t = 4.5 ms for the pure-Ne discharge, and P rad ≈ 1 GW at t = 9 ms for the mixed-D-Ne discharge, where t = 0 is determined by the time of first observed light from the pellet. The TPF is estimated to be reduced as the pure-Ne pellet travels deeper into the plasma from TPF ≈ 1.6 at t = 1.3 ms to TPF ≈ 1. near the end of the thermal quench. The TPF of the mixed-D-Ne pellet is estimated to be TPF ≈ 1.6 for the majority of the thermal quench.
Traces of M3D-C1 simulations of these plume representations are shown in figure 7, with the slower pure-Ne plume in blue and the faster mixed-D-Ne plume in magenta. Both plumes show similar thermal quench times (figure 7(a): ∼t TQ = 3.5 ms for mixed-D-Ne versus t TQ = 3.75 ms for pure Ne) and near-identical onset of a radiation spike at t TQ ∼ 2.7 ms ( figure 7(b)). This is inconsistent with what was observed experimentally, where the mixed-D-Ne pellet had a significantly longer thermal quench. While it is essential to incorporate this effect for a more realistic simulation, adding a pre-pellet injection gas would only accelerate the thermal quench time. Therefore, it is not likely to be the primary cause of the discrepancy between simulations and experiments. One possible explanation for this discrepancy is the MHD simulations with the broad toroidal peaking compared to the experimental cloud width. This broad toroidal peaking reduces the radial shift of the ablated cloud due to the ∇B-drift effect [16,17]. This effect could move the neon away from the hot plasma and delay the thermal quench. This drift has been observed to be small for small neon mixture pellets due to neon radiating cooling the local pressure bump of the plasmoid [18]. However, in these SPI MHD simulations with a large amount of neon, the T e drops to a very low temperature (∼10 eV) behind the pellet [9]. So the plasmoid could drift into a very cold lowpressure region of the plasma. Thus, even with the pressure in the plasmoid rapidly dropping due to radiation, the plasmoid could potentially still be drifting significantly and improve simulation-experimental comparisons.
One sees very different dynamics in the two thermal quenches when the same quantities are plotted versus penetration depth, as can be seen in figures 7(d) and (e). Here we define ρ abl by mapping the average shard location, weighted by each shard's ablation rate, to the normalized toroidal ρ of the initial equilibrium 6 . While the pure-Ne plume only travels to ρ abl =∼ 0.72 before the radiation spike begins, the mixed plume travels ∼0.53. This corresponds to the pure-Ne plume reaching just inside the q = 2 surface before the onset of MHD instability, while the mixed plume penetrates all the way to the q = 5/4 surface before exciting significant MHD activity. As seen in figure 8, while both modes excite multiple toroidal harmonics in the magnetic energy and are dominated by n = 1 activity, there is a proportionally larger n = 2 component in the mixed-pellet MHD crash, indicating greater mode coupling when the plume reaches lower q surfaces. Finally, figures 7(c) and (f ) show the change in the total number of electrons in these simulations, with the mixed pellet increasing the density much more quickly due to increased ablation and ionization. This is due to a combination of the significantly higher ablation rate for deuterium pellets, the increased surface area with more shards, and a higher incident temperature for the mixed plume due to its deeper penetration.
In order to separate the effects of pellet speed and pellet composition, additional M3D-C1 simulations with the speeds of the plumes swapped were performed. These are also plotted in figure 7, with a pure-Ne plume at 300 m s −1 in cyan and a mixed-D-Ne plume at 150 m s −1 in red. At the slower speed, the quench dynamics for the two compositions are nearly identical. Each plume penetrates just inside the q = 2 surface before the onset of macroscopic MHD instability. While the slow, mixed simulation terminated early due to numerical instability from lack of resolution with the current profile becoming grid to grid just before the termination, one can see the beginnings of the radiation spike (figure 7(e)) and MHD instability (figure 8) in that simulation. At the higher speeds, the dynamics are very different between the two compositions. The pure-Ne pellet is unable to penetrate as deeply as the mixed-D-Ne pellet, though still somewhat farther than the slow pellets, approaching the q = 3/2 surface before the onset of MHD instability. This deeper penetration does lead to increased mode coupling between the n = 1 and n = 2, as seen in figure 8. There are likely two causes for the earlier onset of MHD instability with the fast neon pellet compared to the fast mixed pellet. The first is the significantly higher radiation with the fast neon pellet, due to the larger quantity of deposited neon, leading to increased pressure gradients. At the slower speeds, this seems not to matter; the pellet is moving slowly enough inward that even the 5% neon in the mixed pellet is enough to effectively radiate the thermal energy available. The second effect is the higher effective charge number in the plasma with the fast pure-Ne plume compared to the fast mixed-D-Ne plume, which would increase the resistivity and allow for more rapid current diffusion. This could also accelerate the time scale to MHD instability. We can conclude, therefore, that there is a competition of time scales between the rate of radiative dissipation and the movement speed of the pellet. At low speeds, the radiative dissipation dominates and the pellet composition has little effect on quench dynamics. At high speeds, however, pellets with a smaller fraction of radiative material can penetrate much deeper into the plasma, greatly altering the quench dynamics.
NIMROD simulations of the same realistic plumes (i.e. the slow Ne plume and the fast mixed-D-Ne plume) were also performed and found similar predictions in both reductions of stored energy and radiation to that of M3D-C1. The blue and magenta curves in figure 9 show the predicted traces versus time and pellet front for SPI with pure-Ne and mixed-D-Ne pellets, respectively. These curves are equivalent simulations to the M3D-C1 simulation blue and magenta curves of the right column of figure 7. NIMROD predicts a slightly slower thermal quench time/distance with E th reduced by 85% at a front of ρ abl = 0.6 compared to M3D-C1 which had a completed thermal quench at ρ abl = 0.68. Additionally, with the slower thermal quench, the NIMROD simulation has a peak radiation level of 3.5 GW for the pure-Ne pellet, which is slightly lower than M3D-C1 peak radiation of 4.4 GW. The mixed-D-Ne pellet simulations are terminated before the completion of the thermal quench due to a numerical instability. Nevertheless, the mixed-D-Ne simulation before the numerical instability is qualitatively similar to the M3D-C1 result, finding that the pellet travels much deeper into the plasma compared to the pure-Ne pellet, traveling 0.6 m before the thermal energy is reduced to 2.4 MJ.
NIMROD simulations find increasing the amount of ablated material reduces the thermal quench time of both the mixed-D-Ne pellet and the pure-Ne pellet. The purple and orange curves in figure 9 show NIMROD simulations with the ablated material increased by a factor of ten. The number of electrons added to the simulation from the pellet is shown in the bottom panel of figure 9. While the ablation rate was increased by a factor of ten, the amount of free electrons entering the simulation only increased by a factor of 2.07 for the pure-Ne pellet, and by a factor of 5.38 for the mixed-D-Ne pellet when pellet fragments are at ρ abl = 0.75. This is due to the ablated atoms cooling the plasma leading to lower ablation rates and lower fraction of the ablated atoms being ionized. NIMROD predicts a thermal quench t = 3.5 ms for the pure-Ne pellet with artificial 10× material compared to t = 4 ms with 1× material. Along with the faster thermal quench time, the peak radiation increased to P rad = 6.8 GW. The mixed-D-Ne pellet simulation also has a faster thermal quench time with Trace of E th (top (a), (d)), P rad + P ion (middle (b), (e)), and δNe (bottom (c), (f )) for NIMROD realistic SPI plume simulations into JET target discharge 95707. Quantities on the left are plotted versus time and quantities on the right are plotted versus pellet front. Simulation with 1× material using t = 50.544 s equilibria are in blue for pure-Ne and magenta for mixed-D-Ne pellet. Simulation with 10× material using t = 50.875 s equilibria are in purple for pure-Ne and orange for mixed-D-Ne pellet. the increased ablation rate with E th being reduced by 2/3 at t = 2.1 ms up until this simulation is terminated early by a numerical instability.
The TPF = P rad (ϕ = 0)/P rad (n = 0) is predicted to be at levels consistent with experimental observations when the radiated power is above P rad = 0.25 GW, and is reduced as the pellet penetrates deeper into the plasma, as shown in the top panel of figure 10. NIMROD simulations deposits the plume toroidal with a Von Mises distribution with a TPF = 3.63. At the start of the simulations, the simulated TPF is at a level near the deposited SPI plume. This is because the radiating density has not had any time to relax. As the pellet penetrates deeper, the TPF drops to TPF = 2 at t = 1.5 ms, and drops further to TPF = 1.7 at t = 3.5 ms. This reduction in TPF is due to the radiating density having more time to relax. Additionally, when the rate of ablation is increased by a factor of ten, the TPF is decreased shifting it closer to the experimental fit with TPF = 1.6. Like the pure-Ne pellet simulations, the mixed-D-Ne pellet simulations also have a TPF that is at a level near the deposited SPI plume of TPF ≈ 3.65, and a quick reduction of the TPF as time increases. The bottom panel of figure 10 shows the NIMROD simulated TPF versus time for the pure-Ne. At t = 1.25, the simulations predict a similar TPF to the experimentally fitted TPF = 1.6. Interestingly, the radiation becomes nearly symmetric for the artificial 10× material mixed-D-Ne pellet simulation between t = 1.5-2 ms. This is due to the large

Conclusions
Nonlinear 3D MHD simulations using both NIMROD and M3D-C1 show typical SPI disruptions for JET laying the ground work for more-sophisticated validative and predictive modeling of SPI in JET. Switching from a single monolithic to pencil beam to realistic fragment representations modestly affects the thermal quench time. Likewise, a scan in viscosity shows a modest change in the simulations where reducing viscosity increases MHD activity and decreases the thermal quench time by approximately 10%. Realistic cloud of fragments modeling shows that a mixed-D-Ne pellet travels deeper into the core before thermal quench than a pure-Ne pellet. At the slower pellet fragment speed of 150 m s −1 , the pellet is found to be moving slowly enough inward that even the 5% neon in the mixed pellet is enough to effectively radiate the thermal energy available, leading to the pure-Ne and mixed-D-Ne pellet having nearly identical reduction of thermal energy before the onset of MHD activity. Radiation toroidal peaking is reduced as the pellet travels deeper into the plasma, but can have a secondary peak near the end of the thermal quench with the radiation toroidal peaking predicted to be at levels consistent with experimental observations (TBF = 1.6-2.4) when the radiated power is above 1 GW.