Plasmoid drift and first wall heat deposition during ITER H-mode dual-SPIs in JOREK simulations

The heat flux mitigation during the thermal quench (TQ) by the shattered pellet injection (SPI) is one of the major elements of disruption mitigation strategy for ITER. It’s efficiency greatly depends on the SPI and the target plasma parameters, and is ultimately characterised by the heat deposition on to the plasma facing components. To investigate such heat deposition, JOREK simulations of neon-mixed dual-SPIs into ITER baseline H-mode and a ‘degraded H-mode’ with and without good injector synchronization are performed with focus on the first wall heat flux and its energy impact. It is found that low neon fraction SPIs into the baseline H-mode plasmas exhibit strong major radial plasmoid drift as the fragments arrive at the pedestal, accompanied by edge stochasticity. Significant density expulsion and outgoing heat flux occurs as a result, reducing the mitigation efficiency. Such drift motion could be mitigated by injecting higher neon fraction pellets, or by considering the pre-disruption confinement degradation, thus improving the radiation fraction. The radiation heat flux is found to peak in the vicinity of the fragment injection location in the early injection phase, while it relaxes later on due to parallel impurity transport. The overall radiation asymmetry could be significantly mitigated by good synchronization. Time integration of the local heat flux is carried out to provide its energy impact for wall heat damage assessment. For the baseline H-mode case with full pellet injection, melting of the stainless steel armour of the diagnostic port could occur near the injection port, which is acceptable, without any melting of the first wall tungsten tiles. For the degraded H-mode cases with quarter-pellet SPIs, which have 1/4 total volume of a full pellet, the maximum energy impact approaches the tolerable limit of the stainless steel with un-synchronized SPIs, and stays well below such limit for the perfectly synchronized ones.


Introduction
The Thermal Quench (TQ) heat flux mitigation by impurity Shattered Pellet Injection (SPI) is an important part of the ITER Disruption Mitigation System (DMS) strategy, the aim of which is to evenly redistribute the thermal energy stored within the plasma onto the Plasma Facing Components (PFCs), so that localized heat deposition on the divertor or the first wall and the consequential material damage could be avoided during the TQ process [1].The TQ mitigation efficiency thus greatly depends on both the radiated energy fraction, as well as the asymmetry of the radiation heat flux onto the first wall.
There have been intensive studies, both experimentally [2][3][4][5][6] and numerically [7][8][9][10], on both the total radiated energy and the radiation power density within the plasma volume after SPIs.It has been found in general that an apparent toroidal asymmetry in the radiation power density exists for the single SPI cases, especially in the early injection phase.The radiation power density profile tends to relax as the TQ proceeds, and the Toroidal Peaking Factor (TPF) of the radiation power density at the time of the maximum radiation power tends to be smaller compared with the maximum TPF during the whole TQ phase in both experimental observations [3] and numerical simulations [8].
Although the aforementioned studies already provide valuable insight into the asymmetric TQ mitigation process, the ultimate criterion of the disruption mitigation efficiency should be the heat flux onto the first wall and its accumulated energy impact, as the latter is directly related to the wall material temperature rise and thus the wall damage during disruption [11].The recent change of the ITER first wall material to tungsten offers more tolerance to such transient heat flux [12], however, should the accumulated energy impact exceed a certain threshold, even tungsten could experience significant melting [11].Meanwhile, the stainless steel cover for the port window suffers from a lower energy impact threshold which is about 1/3 of the tungsten threshold, thus is more susceptible of melting, although such melting is acceptable as long as the melt threshold is not dramatically exceeded for the stainless steel.[13,14].Furthermore, even without exceeding the tolerable maximum energy impact, transient heat fluxes could already produce structural modification on tungsten surface which could affect the mechanical and thermal properties of the first wall [15][16][17][18][19]. Hence it is desirable to obtain the heat flux and energy impact distribution onto the first wall after SPIs to assess the efficiency of the TQ mitigation and to validate the ITER PFC heat tolerance specifications.
To take a first look into the aforementioned heat deposition behaviour in ITER Hmode plasmas, we carry out 3D non-linear reduced MHD simulations using the JOREK code [20] and its collisional-radiative impurity module [21].Two sets of equilibria are considered, a "baseline" H-mode and a "degraded" one.The latter means to take into account the natural confinement degradation in the disruption precursor phase, which significantly reduces the thermal energy before the triggering of the DMS [22].Two sets of dual-SPI configurations are also considered, one "full pellet" with small neon mixture ratio and one "quarter-pellet" with the same amount of neon atoms but one fourth the hydrogen content.For the full pellet dual-SPI into baseline H-mode case, a strong plasmoid drift along the major radial direction is found to occur due to strong local over-pressure within the newly ablated plasmoid cloud which induces polarization and subsequent E × B drift towards the Low Field Side (LFS) [23].Such drift and its accompanying MHD response would undermine the assimilation when fragments are injected from the LFS mid-plane, as well as cause strong boundary heat flux as the edge becomes stochastic.These detrimental behaviours could be mitigated if the higher neon mixture ratio quarter-pellet is used so that the local over-pressure is suppressed by local radiation [24] or by considering the degraded H-mode.Indeed, the quarter-pellet SPIs into degraded H-mode shows good impurity assimilation and radiation fraction.Moreover, the radiation heat flux onto the first wall is obtained by the Raysect/CHERAB code suite integrated within Integrated Modelling & Analysis Suite (IMAS) [25,26] and their energy impact is calculated.It is found that the maximum energy impact for the full pellet dual-SPI into the baseline H-mode case exceeds the tolerable limit of the stainless steel near the injector, but stays well below that of the tungsten.For the quarter-pellet dual-SPI into the degraded H-mode, the maximum energy impact only exceeds the stainless steel limit slightly even when there is 1ms delay between the injectors, while it stays well below said limit if the synchronization is perfect.
The rest of the paper is arranged as the following: In Section 2, the two target equilibria are introduced, and our basic assumptions, simulation setup as well as the SPI configurations and parameters are presented.We also introduce our method of calculating the accumulated energy impact from a time-varying heat flux distribution onto the first wall.In Section 3, the plasmoid drift behaviour during H-mode SPIs and its mitigation are investigated, and the possible consequences of such drift are discussed.The radiation heat flux and its energy impact onto the first wall for both the baseline and degraded H-mode cases are shown in Section 4, and their implication to the PFC damage is discussed.Finally, conclusion and further discussions on the TQ mitigation efficiency with neon SPIs are presented in Section 5.

Simulation setup and basic assumptions
The 3D non-linear reduced MHD version of JOREK with collisional-radiative impurity model [20,21] is used for this study, and our governing equations are the same as that described in Ref. [9] without modifications, so we shall not repeat them here.The target equilibria are one baseline H-mode chosen from the ITER reference scenarios [27], and one degraded H-mode constructed from the baseline H-mode equilibrium by increasing the perpendicular heat conduction and diffusion artificially thus producing a H-L back transition, mimicking the naturally occurring confinement degradation before the onset of the TQ [22].Both cases have the same toroidal magnetic field at B T = 5.3T and the total plasma current at I p = 15M A. The initial equilibrium profiles, including the electron temperature, density, current density and the safety factor q profiles for both target plasmas are shown in Fig. 1.The initial thermal energy is 190M J for the degraded H-mode case and 370M J for the baseline H-mode case.
The grid we use is presented in Fig. 2, where the boundary approximately represents the position of the thin wall equivalent of the ITER blanket modules [28].We assume the boundary to be ideally conducting in this study.We include 12 toroidal harmonics from n = 0 to n = 11 in our simulations, and our ablative density source is artificially elongated along the toroidal direction due to limitation in the toroidal resolution.This artificial toroidal spread with ∆ϕ N G = 0.3 radian, corresponding to 2.46m to 1.86m as the fragments travel from the edge to the core, enable us to keep the Fourier components up to the first order in strength with the following density deposition shape: Here R f , Z f and ϕ f are the coordinates of the ablating fragment.We choose ∆r N G = 0.08m in this study, but wish to emphasize that the overlapping between the ablation clouds from different fragments tends to create an envelop of plasmoid, so that our result is insensitive to the exact choice of ∆r N G so long as it is no larger than the plume width.
The exact value of the density source is determined by the Neutral Gas Shielding (NGS) ablation model [29], the detailed implementation is described in Ref. [9].The density generated by ablation would then relax by parallel convection driven by the pressure gradient along the field lines, following the continuity equation.The two temperature model described in Ref. [9] is used, where we track the evolution of the electron and the ion temperatures separately.The Braginskii parallel heat conduction [30] is used for each species respectively.We impose an upper limit of this parallel heat conduction corresponding to the value of the Braginskii electron thermal conductivity at T e ≃ 1.6keV , beyond which the heat conductivity would not rise any more.This is an attempt to mimic the flux limit of the free streaming electrons and to prevent the overestimation of Braginskii thermal conduction in the high temperature regime.With the initial density, the perpendicular particle diffusion is set to be constant at 6m 2 /s and the perpendicular heat conduction at 10m 2 /s.The choice of these coefficients are meant approach the gyro-Bohm scaling [31] while still maintain numerical stability.We do not expect the exact value to impact our result drastically since the stochastic transport and the macroscopic convective transport dominate over the turbulent one after the injection for our simulations.The perpendicular viscosity is set to 30m 2 /s and the parallel one to 150m 2 /s, these values are chosen to meet the requirement of numerical instability suppression.The Spitzer resistivity [32] with the effective charge contribution [33] included is used for our simulation.Here, we have neglected the neo-classical resistivity contribution from the trapped electrons as the temperature drop and the density rise after the injection would prevent the electrons to complete a bounce period before collision.This is consistent with the previous finding that the mean-free-path of several hundred eV electrons is comparable with the size of the cold, dense plasmoid around the fragments after ITER SPIs [9].We assume the plasma is transparent to the impurity radiation, hence our result might over-estimate the peak radiation heat flux due to neglecting the plasma opacity.100 1ms Table 1.The injection parameters for the SPI considered in this study.Note that for BH-FP-dt0 there exists an asymmetry between the plumes although they are injected at the same time.
The SPIs are carried out from the outer equatorial ports EQ-08 and EQ-17, which locate at R = 8.4m, Z = 0.685m and opposite toroidal positions, along the major radial direction.The vertex angle of the injection cone is 20 degrees.The fragment size distribution follows that of the Statistical Fragmentation model [34,35].Our two sets of SPI parameters are the "1% neon full pellet SPI" dual-SPIs with 2.5 × 10 22 neon atoms and 1.8 × 10 24 hydrogen atoms for each injector, as well as the "5% neon quarter pellet SPI" dual-SPIs with 2.5 × 10 22 neon atoms and 4.5 × 10 23 hydrogen atoms for each injector.For the full pellet case, the pellet is shattered into 300 fragments, while for the quarter pellet case the fragment number is 100.Their characteristic fragment size is κ −1 p ≃ 1.43mm and κ −1 p ≃ 1.30mm respectively.Both cases use 500m/s reference injection velocity with ±40% velocity spread.We assume the fragment velocity distributes uniformly over the spreading range, although in reality the velocity distribution has a correlation with the fragment size [36].We do not expect this to strongly alter the dynamics we investigate here, since one could usually modify the shattering angle to partly cancel this effect.It should also be noted that, for our baseline H-mode case, we enabled the velocity spread only in the EQ-08 plume and set all fragments in the EQ-17 plume with the same velocity value, albeit with correct direction spread so that there exists an asymmetry in the dual-SPI despite the two plumes being injected at the same time.We nonetheless don't expect this to greatly impact our results, since in reality there would always exist some asymmetry between the plumes even with very good synchronization between the injectors.
We have summarized the cases we investigate in this paper in Table 1, where the target equilibria, the neon and hydrogen injection amount, the fragment number and the delay between the injectors are shown.For all cases, the first fragment of the EQ-08 plume arrives on the LCFS approximately at t = 0.38ms and the last one arrives approximately at 0.89ms.All the fragments from the EQ-17 plume in the BH-FP-dt0 case arrive approximately at t = 0.54ms since there is no velocity spread.The arrival time of the fragments from the EQ-17 plume for the rest of the cases is the same with that of the EQ-08 plume, plus the delay between the injectors shown in Table 1.We will mainly investigate BH-FP-dt0, DH-FP-dt0 & DH-QP-dt0 in Section 3, and BH-FP-dt0, DH-QP-dt0 & DH-QP-dt1 in Section 4.

The energy impact
To estimate the wall material temperature rise in an inertially cooled scenario, it is important to calculate the accumulated heat on the wall surface considering both the incoming heat flux and the conductive cooling of the PFCs.For the simple case of 1D heat conduction concerning a semi-infinite slab and constant heat flux q s , the temperature increase ∆T after ∆t time of heat exposure is [11] ∆T Here b = √ κρc, ρ is the material density, c is its specific heat capacity, κ is its heat conductivity.The temperature rise is proportional, via the material-dependent factor b, to the so-called energy impact ∆Q ≡ q s √ ∆t which itself is material-independent.The maximum tolerable energy impact could then be obtained for each individual material [11,14], thus providing insights into the probability of severe material damage for a given period of heat deposition.
When considering a time varying heat flux which is more relevant during transient TQ process, one could use the fundamental solution of the heat equation to obtain an analytical description of the heat accumulation on the slab surface without numerically solving the 1D heat conduction equation.This is equivalent to considering the heat pulse q s (t) as a series of continuously arriving "heat packets", each of which individually undergoes a 1D random walk towards infinity with step-length determined by the material heat conductivity after their arrival.The total heat distribution after a certain time period ∆t is then the summation of the 1D spatial probability distribution of all the packets.For each such heat packet, the fundamental solution on x ∈ [0, +∞) is Here α is the thermal diffusivity and x is the distance along the 1D direction.Since we are only concerned with the temperature rise at the wall surface, x could be practically set to zero, and we obtain the time decay of each individual heat packet to be proportional to ∆t −1/2 .The energy impact integrating over all the heat packets could then be obtained by the following convolution: Here t 0 is the beginning time of the heat pulse.For a constant q s , one naturally recovers the previous result ∆Q(∆t) = q s √ ∆t.Henceforth we will calculate the energy impact using Eq. ( 4) and estimate the potential heat damage based on it.

Plasmoid drift, accompanied heat flux and its mitigation
For the baseline H-mode BH-FP-dt0 studied here, it is found that the high thermal energy reservoir could result in strong local pressure peak around the vanguard fragments, which induces vertical polarization, which in turn results in E × B drift along the major radial direction [23].In our simulations, the polarization could be implicitly obtained through the evolution of the vorticity equation in Ref. [9]: Here ρ is the mass density, u is the stream function To see this, one realizes that the stream function u is associated with the electric potential, and the second term on the RHS of the velocity definition is essentially the E × B drift.Thus the Poisson bracket {P, R 2 } term in Eq., (5), which includes the ∇B drift contribution, would drive the aforementioned polarization in the presence of the local plasmoid overpressure.Such polarization would then in turn result in the drift motion.The above discussed strong drift motion is found when the fragments arrive on the pedestal of the baseline H-mode.It should be noted that due to the asymmetry in our two plumes as mentioned in Section 2.1, the EQ-08 plume arrives on the pedestal first and triggers the strong drift.The resulting drift motion expels significant amounts of injected materials out of the Last Closed Flux Surface (LCFS), undermining the injection assimilation, as can be seen in Fig. 3, which shows to the same toroidal location as the EQ-08 plume.
The local electron over-pressure could be readily identified at t = 0.41ms in Fig. 3(a) as the fragments penetrate the equilibrium LCFS represented by the white contour.This over-pressure region essentially is the envelop of several plasmoid produced by individual vanguard fragments as they enter the pedestal.This high pressure plasmoid envelop is then seen convected outward in Fig. 3(b) due to the major radial drift and finally relaxes outside of the separatrix in Fig. 3(c).Consequentially, the injection assimilation within the Last Close Flux Surface (LCFS) suffers in the early injection phase as is shown in Fig. 4, where the comparison between the mid-plane cut of the n = 0 electron density profile for BH-FP-dt0 and DH-QP-dt0 at various times is shown.The LCFS on the low field side approximately locates at R = 8.2m for both cases.In Fig. 4(a), comparing the black, blue and green lines, it can be seen that the density peak initially occurs near the LCFS, then goes outward into the scrape-off layer.This drift motion, coupled with the associated outgoing heat flux, causes significant density deposition in the open field line region, which is then lost over time.As time goes on however, the fragments ultimately penetrate deeper into the plasma and result in core density rise in the late injection phase, as shown by the cyan and red lines.Note that the density rise in the axis region is due to inward MHD transport, as at that time the majority of fragments has not reached the axis yet.This behaviour is compared with DH-QP-dt0 shown in Fig. 4(b), where no apparent density expulsion is seen and the density peak gradually moves inward.
Such drift motion is accompanied by strong field line stochasticity in the edge region, resulting in a stochastic parallel heat flux, as well as perpendicular convective   heat flux due to the drift transport of the high pressure plasmoid.Poincaré plots of the magnetic field lines at t = 0.41ms and t = 0.46ms for BH-FP-dt0, as well as the simulation boundary and a surface 10cm inside of that boundary which corresponds to the approximate position of the first wall, are shown in Fig. 5, where the field lines that remain well confined after 2000 toroidal turns are marked by red, while those that are lost are coloured according to the distance they travel before they hit the simulation boundary.The figure corresponds to the same toroidal location as the EQ-08 plume.The simulation boundary is shown in purple, and the first wall approximation is shown in green.In Fig. 5(a), the edge region begins to become stochastic, and the "connection length" of the field lines initialized in the edge region to the wall decreases.Another notable feature is the outward deformation of the closed flux surfaces close to the plasmoid position, likely a result of the outward plasma motion.In Fig. 5(b), already significant region around the previous pedestal has become stochastic with short connection length to the wall, and the H-mode pedestal is collapsed.Comparing Fig. 4(a) and Fig. 5(b), it can be seen at t = 0.46ms there is no significant density assimilation within the pedestal, while strong edge stochasticity and thus thermal energy loss already occurs.
As a result of the above-mentioned stochastic and convective transport, the conductive and convective heat flux onto the first wall could be significant.In Fig. 6, the normal component of the combined heat flux normal to the first wall approximation represented by the green line in Fig. 5 is shown for times t = 0.42ms and t = 0.43ms.Comparing Fig. 3 and Fig. 6, it can be seen that as the over-pressured plasmoid drifts towards the wall at t = 0.42ms, a strong helical boundary heat flux structure begins to appear at the EQ-08 injection location centred around ϕ = −2.57.It successively relaxes along field lines and the peak heat flux locations move away from the plasmoid drift location at t = 0.43ms.During this transient heat pulse the perpendicular heat flux dominates over the parallel one.
Both the aforementioned density expulsion and the accompanying transport event are undesirable for TQ mitigation.The H-mode pedestal collapses after the transport event is triggered and over 100M J of thermal energy is prematurely lost before the fragments penetrate deep into the core, compared with the 370M J initial thermal energy.However, the plasma usually would already experience significant confinement degradation due to the growth of precursor modes before the disruption onset [22] and the triggering of the SPIs.Hence it is interesting to see if the reduced thermal energy reservoir of such a "degraded H-mode" could still trigger such strong drift motion, and whether an increased neon mixture ratio could help to suppress the plasmoid drift by radiatively reducing the local over-pressure, and thus mitigate the accompanying outgoing heat flux.Indeed, if we compare the thermal energy loss and the radiated power of BH-FP-dt0 and DH-QP-dt0 as is shown in Fig. 7, DH-QP-dt0 enjoys much better radiation fraction behaviour due to the mitigation of this drift motion.Here the total loss of the thermal energy is represented by the black solid lines, and the total radiated energy is represented by the red solid lines.For BH-FP-dt0, the sudden rise of the discrepancy between the thermal energy loss and the radiated energy beginning at t = 0.45ms is due to the aforementioned heat flux accompanying the drift motion.This heat flux is absent in DH-QP-dt0 as will be shown below.For both cases, the ohmic heating is negligible until the end of the simulation when it begins to pick up due to the thermal energy reservoir getting depleted.The perturbed magnetic and kinetic energies with toroidal harmonics n = 1 to n = 6 for DH-FP-dt0 and DH-QP-dt0 are shown in Fig. 8, where the perturbed magnetic energies are shown in solid lines and the kinetic ones are shown in dashed lines.We identified two MHD peaks in both cases, which are marked by the vertical chained lines in the figure.For both cases, the first peak corresponds to the fragment arrival on the plasma which results in edge stochastization, while the second peak corresponds to the core temperature relaxation.These strong MHD responses correspond to strongest outgoing heat flux and plasmoid drift motion, thus we would investigate the boundary heat flux at these times.
The log 10 of the normal component of the boundary heat flux at the two times when the MHD amplitude reaches its peak values for both cases are shown in Fig. 9. Fig. 9(a) and (b) correspond to two MHD peaks of DH-FP-dt0 and Fig. 9(c) and (d) correspond to that of DH-QP-dt0.For DH-FP-dt0, a small asymmetry in the MHD response initially results in a weak drift motion and mild out going heat flux from the EQ-08 plume position at ϕ = −2.57and t = 0.74ms as is shown in Fig. 9(a).Later on, at the time of the core relaxation, strong heat flux coming out from the X point of the 1/1 magnetic island results in a local high pressure plasmoid close to the EQ-17 plume location at ϕ = 0.57 and t = 1.17ms, which is accompanied by a strong outgoing heat flux shown in Fig. 9(b).The heat flux associated with the second drift is comparable in amplitude with the heat fluxes in BH-FP-dt0 as is shown in Fig. 6.With higher neon mixture ratio, DH-QP-dt0 found negligible plasmoid drift, and almost no accompanying outgoing heat flux even during the MHD peak, as is shown in Fig. 9(c) &  (d).This is responsible for the difference in the radiation fraction between BH-FP-dt0 and DH-QP-dt0 shown in Fig. 7.
To summarize, the arrival of the vanguard fragments on the baseline H-mode pedestal would trigger strong outward plasmoid drift which is accompanied by edge stochastization, resulting in significant density expulsion and a large heat flux onto the wall.This is detrimental to the TQ mitigation.Considering the confinement degradation in the precursor phase as well as using higher neon concentration pellets can effectively mitigate such drift motion and the associated heat flux, thus increasing the radiation fraction.

Radiation heat flux and its energy impact after neon dual-SPIs
Another important issue that has to be addressed is the asymmetry of the radiation heat deposition onto the ITER first wall.Concerns have been raised that if there exists strong asymmetry in the volumetric radiation power density, local parts of the first wall could still receive significant amount of heat flux leading to surface melting of PFCs.For the design of the ITER in-vessel components such as diagnostics or heating systems, the radiation peaking and the resulting heat fluxes to these components have been specified.However, validation of these specifications is still outstanding and the simulations presented here contribute to this validation process. is calculated by the method detailed in Section 2.2, which then is compared to the maximum tolerable energy impact limit to assess the potential damage to the PFCs.This serves as a first demonstration of the JOREK capability to provide input for future more detailed PFC heat damage investigations.It should be noted that directly using the high time resolution JOREK data to calculate the energy impact is impractical due to the enormous data size, so we used lower time resolution data with linear interpolation between time slices to carry out the time integration needed in Eq. ( 4).We mainly focus on three cases in this section: BH-FP-dt0, DH-QP-dt0 & DH-QP-dt1 as shown in Table 1.The overall assimilated neon, assimilated hydrogen within the LCFS, radiated energy fraction f rad , approximate TQ timescale defined as the timescale of the thermal energy falling from 90% to 20% of the initial thermal energy, as well as the maximum local energy impact for those three cases are summarized in Table 2.We are interested in DH-QP-dt0 & DH-QP-dt1 since the quarter pellet dual-SPIs into the degraded H-mode shows good radiation fraction, as is demonstrated in Fig. 7(b).We are interested in the comparison between the perfectly synchronized DH-QP-dt0 and the 1ms delayed DH-QP-dt1, since the mismatch between the arrival time of the fragment plumes means stronger asymmetry in the radiation power density and probably a larger maximum heat flux, and it is of interest to see if the 1ms delayed case would see local radiation energy impact exceeding the tolerable limit of the PFCs.The baseline H-mode BH-FP-dt0 is included to demonstrate the expected strongest radiation heat deposition onto the first wall due to the large radiation power, despite its smaller radiated fraction.
The correlation between the radiation power, the neon ablation rate of both plumes, as well as the n = 1 and n = 2 magnetic energy perturbation for all three cases are shown in Fig. 10, with red lines representing the radiation power, the blue solid lines representing the neon ablation rate of the EQ-08 plume, the blue chained lines representing that of the EQ-17 plume, the black solid lines representing the n = 1 perturbed magnetic energy and the black chained lines representing the n = 2 perturbed magnetic energy respectively.The approximate arrival time of the vanguard fragments on the equilibrium LCFS is represented by the green vertical dashed line for the EQ-08 plume and the magenta vertical dashed line for the EQ-17 plume.For BH-FP-dt0 shown in Fig. 10(a), the very sharp peak in ablation and radiation at approximately t = 0.42ms coincide with the strong plasmoid drift discussed in Section 3, the following peak in the MHD amplitude corresponds to the edge stochasticity shown in Fig. 5.The later radiation peaks coincide with the core temperature relaxation and cooling,  which corresponds a weaker MHD peak just before t = 1.5ms.For DH-QP-dt0 shown in Fig. 10(b), there are two radiation peaks in the early phase of injection at approximately t = 1.1ms and t = 2ms, the later peak corresponds to the core temperature relaxation.These peaks are closely correlated with peaks in ablation and MHD amplitude.This correlation is the result of very nonlinear interaction between the fragments and the MHD mode.The ablation creates locally peaked lowly charged impurity density profile, thus result in strong local radiation cooling.This local cooling and local pressure peak incurs helical current perturbation which destabilize resonant MHD modes.The growing MHD modes result in enhanced outgoing heat flux which further strengthen the ablation in the plume.Thus the strong correlation between the radiation, ablation and MHD amplitude peaks especially in the early injection phase.DH-QP-dt1 shows similar correlation as is shown in Fig. 10(c).
For each case, a few time slices are chosen to show both the instantaneous heat flux and the accumulated energy impact onto the first wall.The time around the radiation peaks shown in 10 are especially of interest, so we choose the time slices to be at the radiation peak for the heat flux analysis and right after the peak for the energy impact analysis for the early injection phase.For completion, a few time slices are also chosen for the late injection phase when the radiation power density distribution is already relaxed within the plasma.The ITER first wall tiles as well as the diagnostic window armour plates are shown in Fig. 11  We first take a look at DH-QP-dt1.The instantaneous heat flux onto the first wall is shown in Fig. 12.Since the EQ-08 plume arrives earlier in this case, the radiation heat flux concentrates around the EQ-08 port in the early injection phase as can be seen from Fig. 12(a) and (b), at t = 0.52ms and t = 0.90ms respectively.These two times correspond to the first two peaks in the radiation power shown in Fig. 10(c).The toroidal elongation of the heat spots in Fig. 12(a) and (b) is due to the toroidally elongated density source shape as we have discussed in Section 2.1.It can be seen that initially the radiation heat flux is almost symmetric around the EQ-08 port in Fig. 12(a), but shifts to the right in Fig. 12(b).This is a result of the interplay between the impurity density distribution and the outgoing heat flux which moves the radiation power density peak away from the fragment location that has been reported previously [8].Ultimately, the radiation power density relaxes over the plasma volume both due to the impurity density relaxation and the arrival of the second plume.As a result, the first wall radiation heat flux becomes more or less uniform in the late phase of the SPI, as can be seen in Fig. 12(c 17 at t = 5.74ms.At this time, despite the high total radiation power as is shown in Fig. 10(c), the maximum radiation heat flux onto the first wall is one order of magnitude lower compared with the early injection phase shown in Fig. 12(a) and (b).Another notable feature is that the heat flux on the armour plate is smaller compared with that of the first wall tiles, which is due to a geometric effect.The energy impact factor at several chosen times for the same case is shown in Fig. 13.Fig. 13(a) and (b) show the energy impact just after the radiation spike at t = 0.56ms and t = 0.94ms, corresponding to Fig. 12(a) and (b), and indeed one can see the energy impact distribution correlates well with that of the transient heat flux.At the end of the second radiation spike at t = 0.94ms, part of the stainless steel armour plate shows the energy impact approaching the maximum tolerable limit of 13M Js −1/2 /m 2 [14].Another feature is that the tiles beside the EQ-08 plate show higher energy impact, consistent with the heat flux behaviour shown in 12(b).Nevertheless, the energy impact on those tiles is well below the maximum tolerable limit for the tungsten at 38M Js −1/2 /m 2 [11].In the late injection phase at t = 5.74ms, despite the high total radiation power, the energy impact onto the first wall is much more uniform and stays well below the stainless steel limit, as is shown in 13(c) and (d) which look at the EQ-08 and EQ-17 port respectively.As a comparison, the first wall energy impact in the early and late injection phase for DH-QP-dt0 is shown in Fig. 14.The energy impact at t = 2.06ms, after the second radiation peak shown in Fig. 10(b), around the EQ-08 and EQ-17 ports is shown in Fig. 14(a) and (b) respectively, while the late time energy impact at t = 5.64ms is shown in Fig. 14(c) and (d) for both locations.Due to the good synchronization between the injectors, both ports see similar energy impact accumulation in the early injection phase, and the maximum energy impact stays well below that of the stainless steel limit, let alone the tungsten one.In the late injection phase, due to the relaxation of the radiation power density, the energy impact becomes rather uniform across the torus apart from some poloidal asymmetry due to the shape of the first wall cross-section, as can be seen in Fig. 14(c) and (d).Again, despite the higher radiation power, the maximum energy impact stays well below the tolerable limit of the stainless steel.
Finally, we take a look at BH-FP-dt0 to demonstrate the highest radiation energy impact in our simulations.This high energy impact is because of both the high total radiation power and the asymmetry in the fragment plume we explained previously.Due to this asymmetry, we expect the radiation to concentrate around the EQ-08 port since the vanguard of its plume would first reach the plasma and trigger the plasmoid drift motion discussed in Section 3. Indeed, as is shown in Fig. 15(a) and (b) at t = 0.43ms, after the sharp radiation peak caused by the vanguard fragments arriving on the baseline H-mode pedestal, significant energy impact is accumulated around the EQ-08 port while the energy impact around the EQ-17 port is negligible.
The energy impact of this radiation pulse already exceeds the tolerable limit of the stainless steel, thus we expect to see melting of the EQ-08 armour plate even without considering the conductive/convective heat flux accompanying the plasmoid drift.The maximum radiation energy impact is still well below the tungsten limit, however.As time goes on and the fragments fly further inward, we see the radiation asymmetry turns to the other side of the torus around the EQ-17 port, as is shown in Fig. 15(c) and (d) at t = 2.00ms.The radiation energy impact also concentrates around the injection port at this time, but shows stronger helical feature, representing the parallel relaxation of the impurity density.Again, the maximum energy impact exceeds the stainless steel limit, but stays well below the tungsten one.

Conclusion and discussion
JOREK simulations with collisional-radiative impurities are carried out for dual-SPIs into ITER H-mode plasma to investigate the heat deposition during disruption mitigations.Two sets of equilibria are considered, one baseline H-mode following the default ITER H-mode scenario, and one degraded H-mode where an artificial H-L back transition is induced to mimic the confinement degradation in the disruption precursor phase.Into these two equilibria, two sets of dual-SPI configuration are simulated, the full pellet case with small neon fraction and the quarter pellet case where the neon amount is the same but the hydrogen amount is reduced to one quarter of the full pellet case.
It is found that with baseline H-mode and full pellet, a large local plasmoid over-pressure occurs near the vanguard fragments, which induces polarization and consequential drift motion outward along the major radius.This drift motion causes evident density expulsion and is accompanied by edge stochasticity and strong outward heat flux, both detrimental to the disruption mitigation.The transient heat flux is so large that even with expose time on the order of 1µs, the resulting energy impact could still be on the order of 10 8 M Js −1/2 /m 2 , exceeding the tolerable limit of the tungsten.By considering the pre-disruption confinement degradation and using the higher neon fraction quarter pellet, such local over-pressure could be reduced, thus mitigating the drift motion and the accompanying heat flux.In these scenarios, the transient heat flux is significantly diminished and the radiation fraction during the TQ is improved.
The radiation heat flux and its accumulated energy impact onto the first wall is obtained by ray-tracing from the volumetric radiation power density distribution provided by JOREK.Three cases are investigated: the baseline H-mode with full pellet case and two degraded H-mode with quarter pellet cases with perfectly synchronization and 1ms delay between the injectors respectively.It is found that the radiation heat deposition tends to concentrate around the injecting port in the early injection phase, while it spreads over larger areas towards the end of the TQ.With the quarter pellet, degraded H-mode and 1ms injection delay, the maximum radiation energy impact approaches the tolerable limit of the stainless steel plate while staying far below the tungsten limit.Good synchronization is found to further relax the energy impact distribution and reduce its maximum value, such that it never reaches the stainless steel limit for the perfect synchronized degraded H-mode case.These results are reassuring since they suggest the degraded H-mode with quarter pellet case could achieve high radiation fraction without the risk of melting the first wall.In the baseline H-mode case, despite the strong radiation power, its energy impact still stays well below the tungsten limit, although it exceeds that of the stainless steel substantially, hence one would expect melting of the diagnostic window armour plates.However, even if such melting occurs, it would not impede continued operation of ITER, as it has been found in previous studies with 22M Js −1/2 /m 2 heat pulses that the transient surface melting would result in a surface roughening increasing at a rate of 1-2 micrometers per pulse, but no loss of material [13].Another interesting observation is that the maximum energy impact is approaching the melting limit of the beryllium at 28M Js −1/2 /m 2 [14], so that the new ITER tungsten tiles provide more margin against radiation heat deposition compared with previous beryllium ones.
A few assumptions in our investigations might impact our results presented above.First of all is our toroidally elongated ablation density source, which result in enhanced impurity density spreading along the toroidal direction.One could expect that a more realistic, and expensive, ∆ϕ N G would induce stronger impurity density peak and radiation power density peak, thus more localized radiation heat flux.This would especially be important if the radiating source is very close to the wall as in the early injection phase, or if the plasmoid drift causes the radiating plasmoid to move outward towards the wall.Another consequence of this elongated toroidal density source is that it tends to result in less concentrated pressure peak since the density is more spread out, thus reducing the drive for the plasmoid drift.It is however hard to say if a more realistic toroidal elongation would result in stronger drift, as the accompanying stronger radiation would tend to reduce the local pressure, thus counter-act the pressure peaking effect discussed above.This is left for future studies.Furthermore, we also did not consider the opacity of the plasma.The high density plasmoid around the fragment plume could absorb significant portion of the radiation and re-emit them at a later time [37], reducing the radiation peak.Resolving these two competing factors is part of our future plan.Last, our choice of viscosities might impact the TQ duration, as is demonstrated in Ref [10].The current values we chose are the minimum value we could use without rampant numerical instabilities.Using smaller viscosity might results in more rapid MHD growth and quicker TQ.However, as is shown in Ref [10], changing the viscosity from 2000m 2 /s to 500m 2 /s only reduces the TQ duration from 2.9ms to 2.6ms in their viscosity scan.Hence we don't expect drastic changes to our conclusion should a more realistic viscosity is used.
Nevertheless, our study reported above serves as a first look into the radiation heat deposition on the ITER first wall during its H-mode TQ mitigation process using nonlinear 3D simulation, which is part of ongoing validation work for the ITER PFC heat load specification.This study also provides input for more detailed heat transport and material damage studies in the future.

Figure 1 .Figure 2 .
Figure 1.The initial equilibria of (a) the baseline H-mode and (b) the degraded H-mode, ψ is the normalized poloidal flux.

Figure 3 .
Figure 3.The electron pressure profile at the same toroidal location as the EQ-08 plume for BH-FP-dt0 at (a) t = 0.41ms, (b) t = 0.42ms and (c) t = 0.46ms.The white contour indicates the position of the equilibrium LCFS.The over-pressure as a result of the vanguard fragments arriving on the pedestal as well as its fast outward drift and expulsion can be identified.

Figure 4 .
Figure 4.The mid-plane cut of the n = 0 electron density for (a) BH-FP-dt0 and (b) DH-QP-dt0 at various times.The high pedestal temperature and the consequential plasmoid drift in BH-FP-dt0 results in significant density deposition outside of the pedestal in the early injection phase as is shown by the black, blue and green lines in (a).The LCFS on the low field side approximately locates at R = 8.2m.

Figure 5 .
Figure 5. Poincaré plots of the magnetic field lines at time (a) t = 0.41ms and (b) t = 0.46ms for BH-FP-dt0.The field lines remain well confined after 2000 toroidal turns are marked by red, while those that are lost are coloured according to their "loss length".The outward deformation of the flux surfaces on the mid-plane in (a) is likely due to the outward plasma motion.

Figure 6 .Figure 7 .
Figure 6.The normal component of the combined conductive and convective heat flux for BH-FP-dt0 at time (a) t = 0.42ms and (b) t = 0.43ms across a surface 10cm within the simulation boundary, which represents the approximate location of the thin wall equivalent of the blanket module.The horizontal axis is the toroidal angle and the vertical one is the poloidal angle.The EQ-08 port locates approximately at ϕ = −2.57.During this transient heat pulse the perpendicular heat flux dominates over the parallel one.

Figure 8 .
Figure 8.The MHD spectrum of (a) DH-FP-dt0 and (b) DH-QP-dt0.The solid lines represent the perturbed magnetic energy and the dashed lines represent that of the kinetic energy.The red chained lines corresponds to the MHD peaks at the time of which we would investigate the boundary heat flux distribution.

Figure 9 .
Figure 9.The normal component of the combined conductive and convective heat flux of DH-FP-dt0 at times (a) t = 0.74ms and (b) t = 1.17ms across a surface 10cm within the simulation boundary, as well as that of DH-QP-dt0 at time (c) t = 1.09ms and (d) t = 1.90ms.All four time slices coincide with peaks of MHD amplitude in their respective cases.The EQ-08 port locates approximately at ϕ = −2.57and the EQ-17 port locates at ϕ = 0.57.

Figure 10 .
Figure 10.The correlation between the radiation power, the neon ablation rate of both plumes, as well as the n = 1 and n = 2 magnetic energy perturbation for (a) BH-FP-dt0, (b)DH-QP-dt0 and (c)DH-QP-dt1.The approximate arrival time of the vanguard fragments on the equilibrium LCFS is represented by the green vertical dashed line for the EQ-08 plume and the magenta vertical dashed line for the EQ-17 plume.For the DH-QP-dt0 case, the two arrival time overlap with each other.

Figure 11 .
Figure 11.The ITER first wall (white) and the diagnostic port armour (red).The EQ-08 port locates at X = −7.05mand Y = 4.57m, it is the third port from right visible in the figure, marked by the green rectangle.The EQ-17 port locates at X = 7.45m and Y = −3.84m, on the other side of the torus.
, where the tiles are shown in white and the armour plates are shown in red.The EQ-08 port locates at X = −7.05mand Y = 4.57m, it is the third port from right visible in the figure.The EQ-17 port locates at X = 7.45m and Y = −3.84m, on the other side of the torus.

Figure 12 .
Figure 12.The first wall heat flux of DH-QP-dt1 at (a) t = 0.52ms looking at port EQ-08, (b) t = 0.90ms looking at port EQ-08, (c) t = 5.74ms looking at port EQ-08 and (d) t = 5.74ms looking at port EQ-17.By the time of t = 5.74ms, the radiation heat flux already becomes relaxed, and despite the high total radiation power, the radiation heat flux is low.
under the implementing agreement IO/IA/19/4300002053.This work was supported in part by the Swiss National Science Foundation.Part of this work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 -EUROfusion).The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI).Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI.Neither the European Union nor the European Commission nor SERI can be held responsible for them.This work is carried out partly on Tianhe-3 operated by NSCC-TJ and partly on the ITER cluster.

Table 2 .
The neon and hydrogen assimilation, radiated fraction, TQ time and maximum local energy impact for three cases considered in this section.
To investigate such radiation heat deposition, we utilize the IMAS integrated RaySect/CHERAB code suite to obtain the wall heat flux via ray-tracing the radiation power density obtained by JOREK onto the realistic ITER first wall tiles and the stainless steel armour plate of the diagnostic port windows.Once the heat flux distribution is obtained, the accumulated energy impact