The mechanism of the global vertical force reduction in disruptions mitigated by massive material injection

Disruptions lead to a rapid loss of thermal and magnetic energy and can cause large heat loads, mechanical forces


Introduction
Tokamak experiments regularly observe a fast loss of the thermal and magnetic energy due to disruptions [1][2][3].Large heat loads and electro-magnetic (EM) forces result from the loss of the stored energy due to losses on stochastic field lines, Ohmic heating and induction of current in the surrounding conductive structures.In high current devices like JET, ITER or DEMO, these events have to be avoided or mitigated to avoid life time reduction of components [4].
The issue of vertical stability is closely linked to this topic.In one type of disruption, a loss of the vertical plasma position control leads to a so called vertical displacement event (VDE) [5].In hot VDEs, the contact with the plasma facing components and the associated change of the edge safety factor [6] trigger a thermal quench (TQ), resulting in a current quench (CQ) phase due to the enhanced resistivity.These events are accompanied by currents in the scrape-off layer (SOL) that arise due to the conservation of poloidal and toroidal flux [7], which close through the conductive structures around the plasma and are named halo currents.Together with eddy currents that are induced in conducting structures by time variations of the magnetic field, they contribute to the vertical force on the machine [5,8].In case of toroidal asymmetries, also large horizontal forces can arise [9] which cause a large sideways displacement in JET [10].
VDEs are also observed after major disruptions that result from various causes including mode locking or density limit violation.In such cases, vertical position control is lost as a result of the current decay [5], thus, they are referred to as cold VDEs due to the loss of thermal energy before the vertical movement starts.If the current loss rate is faster than the wall time, the plasma axis position becomes a function of the plasma current [11,12].In addition, a vertical instability occurs when I p reduces below a threshold, where wall currents cannot provide sufficient stabilization below [5,13].
Mitigation techniques aim to reduce the heat loads and EM-forces without facilitating the generation of runaway electrons (REs) [14].The most promising mitigation technique relies on the injection of large quantities of hydrogen to dilute the plasma and prevent the creation of REs and impurities to radiate a large fraction of the stored thermal and magnetic energy in the plasma.Common methods are the injection of a large amount of material as gas (massive gas injection, MGI) or in the form of pellets.The ITER disruption mitigation system (DMS) relies on shattered pellet injection (SPI) from several injection locations [15,16].SPI is experimentally tested in the tokamaks DIII-D [17], HL-2A [18], KSTAR [19], JET [20], J-TEXT [21] and ASDEX-Upgrade [22] at present.Experiments with MGI in JET [23,24], DIII-D [25], Alcator C-mod [26], KSTAR [27] and ASDEX-Upgrade [28] have been routinely carried out to assess the mitigation efficiency against heat loads as well as EM forces and to protect the machine against the consequences of disruptions.The SPI technology is planned for ITER as it promises a higher mass deposition [29] than MGI and mass delivery closer to the core.Recently, an SPI system has been installed in ASDEX Upgrade [22] and the radiation diagnostics by fast cameras and additional bolometers have been improved in support of ITER studies [30].These experiments assess the impurity assimilation, radiation fraction, asymmetries and mitigation efficiencies.Observations over various devices show a reduction of the heat loads and EM-forces compared to non-mitigated disruptions for both SPI and MGI.However, fundamental understanding regarding the mechanisms leading to a reduction of the vertical force is still limited and reliable predictions of the amount by which the forces are mitigated is missing.
In this article, we present a theory for the reduction of vertical forces by massive material injection, carry out axisymmetric simulations with the non-linear extended MHD code JOREK to test our hypothesis, validate the assumptions by comparing to experimental data from ASDEX Upgrade (AUG) and JET, and perform predictions for future ITER experiments.The rest of the article is structured as follows.First, the theory behind the force reduction in mitigated disruption is laid out in section 2, before describing SPI mitigated VDE experiments in AUG in section 3. Third, the computational model of the MHD code JOREK is described and its results compared with the experimental data in section 4.This model is then also applied to ITER cases in section 4.5.Lastly, the findings are summarized and discussed in section 5.

Theory
In this section, we recapitulate the origin of the vertical wall force and discuss the theoretical limit for its reduction.At the end of the section, we propose a mechanism to explain the observed reduction of the current moment, which is based on the stabilization provided by halo currents.

Vertical wall force and vertical current moment.
As it was shown by Miyamoto [31], the vertical force on the vessel F z can be expressed as: where F p,c and F v,c are the forces produced by the poloidal field (PF) coils on the plasma and on the vacuum vessel (VV), respectively.The first term is the source of the wall force and can be expressed as F p,c ∝ I p ∆Z curr [5].The second is a damping term that describes the shielding from the wall currents on the time scales of the resistive decay in the VV.When the total time of the disruption is comparable to the wall time τ w or smaller, this effect has to be taken into account.The expression I p ∆Z curr is similar to the change in the vertical current moment M IZ .The quantity M IZ is defined by: where Z is the vertical coordinate and j ϕ is the plasma toroidal current density (including the halo region).The vertical force is therefore closely related to ∆M IZ and even identical if Z = 0 in the initial equilibrium.In the remainder of the manuscript, we refer to I p ∆Z curr as the change of the vertical current moment.
Theoretical limit to the mitigation efficiency.The mitigation techniques discussed above aim to transfer the majority of the thermal and magnetic energy to radiation by injection of impurities.The vertical wall force is reduced accordingly as the plasma displaces at lower I p , reducing the change in M IZ .However, according to theory, a vertical instability develops at a large external field curvature or when the plasma current decreases below a certain limit [5,13].Moreover, the displacement of the magnetic axis is a function of the plasma current for fast CQs (τ CQ < τ w ) [11,12] as it is the case for mitigated disruptions in ITER.Therefore, the plasma comes into contact with the wall at the same value of I p (in the absence of halo currents) that is given by the threshold I p,th .Consequently, it should not be possible to reduce ∆M IZ below a value given by I p,th : where ∆z max is the vertical distance between the wall and the equilibrium plasma position and I p,th is the plasma current below which the plasma becomes vertically unstable.Due to the relation between F z and ∆M IZ (1), this should also limit the mitigation efficiency for the global vertical force.However, this contradicts experimental observation that mitigation can greatly reduce the vertical force on the structures [23,25,28].Furthermore, disruption mitigation experiments with SPI and MGI commonly observe that the current centroid, i.e. the geometric center of the current distribution, becomes stationary after the TQ and moreover, the vertical motion even seems to reverse its direction in particular cases (see [32][33] and [34]).

Magnetic axis and current centroid.
To solve this apparent contradiction, we propose to distinguish between the vertical motion of the magnetic axis (where B p = 0) and the vertical motion of the current centroid.The latter is defined as: where j ϕ is the toroidal current density and the integration is performed over the full plasma region (including the SOL).
While the magnetic axis is related to the center of the closed magnetic surfaces, the current centroid indicates the geometric center of the full plasma current, including the halo region.Figure 1(a) shows the results of a JOREK simulation (see section 4) of a mitigated disruption during a VDE, where the movement of the magnetic axis is accelerated as the plasma current crosses the I p threshold, while the current centroid stays close to the midplane.The mechanism behind this is explained in the following.

Mechanism of current centroid stagnation during the CQ.
When a radiative collapse leads to the TQ, the balance of Ohmic heating and radiation determines the post-TQ temperature [23].In this phase, the temperature of the core can become comparable to the SOL, and thus, core and SOL resistivities become similar.As a result, halo currents are induced efficiently in this region and the current profile becomes flat, extending beyond the last closed flux surface (LCFS).The diffusion of the currents beyond the LCFS in radiation dominated post-TQ plasmas has already been described in [33] based on observations in DIII-D.This allows the core plasma to move vertically, while the current centroid stays in the midplane as a large fraction of the total plasma current is outside the LCFS.
An example of this behavior is shown in figure 1, showing the time traces (a) of I p and the toroidal halo current as well as the magnetic axis and current centroid and their trajectories (b) in a simulation.At the beginning the current centroid and the magnetic axis displace at the same rate.After the injection, the trajectories diverge and the center of the current stays closer to the midplane, while the center of the plasma core still moves downward.After the LCFS is lost (46.8 ms), the magnetic axis is no longer defined and the current centroid motion reverses due to current diffusion.
As the vertical force depends on the vertical current moment and not on the magnetic axis, its value is determined by the current centroid movement.This allows to mitigate the vertical forces of VDEs, where the maximum of the vertical moment is determined by the displacement at the time of injection.Whereas the magnetic axis can be displaced further as according to theory [11] without influencing F z .In the case of a fast CQ (τ CQ ≪ τ w ), the edge safety factor increases due to the reduction of I p .As a result, the poloidal halo currents can become small despite a large magnitude of the toroidal component due to their relation by q in the halo region q h [7], Another consequence of the increased value of the edge safety factor q 95 is the absence of sideways forces, as they are for example observed in JET [10,35].They are commonly linked to the occurrence of a 1/1 kink mode, which requires a low value of the edge safety factor [36], like commonly observed in hot VDEs.
If the current centroid is stationary during the disruption, one might think it is possible to eliminate the vertical force on the vessel completely.However, due to vessel asymmetries and top-down asymmetries in the external magnetic field, the minimum force is given by the current transferred to the conductive structures during the CQ and their interaction with the vessel.For ITER disruptions, this force is in the order of 10 MN [37].

ASDEX Upgrade experiments
Setup.Recently, an SPI injection system as well as an upgrade on the bolometer and fast camera system has been installed in ASDEX Upgrade, with the goal to perform detailed studies in prospect for ITER [22].The results presented here consist of three SPI mitigated VDEs that were performed with the SPI system to study the mitigation of the loads on the machine for displaced plasmas.The discharges are based on the hot VDE discharge # 39655, that has been described in detail in [38].In order to get a reproducible motion, the VDE was triggered by a downward displacement of 2 cm within 20 ms by the plasma by the control coils.Afterwards, the PF coil currents are frozen, to prevent further influence by the control system.In the discharges # 40955-40957, the initial axis movement coincides almost perfectly.
As the movement of the axis is known from # 39655, the trigger time of the SPI was pre-programmed to initiate a TQ at different phases of the displacement.The pellet consisted of 100% Neon, with a length of 8.5 mm, a nominal diameter of 8 mm and a speed of 242 m s −1 -257 m s −1 .The rectangular 25 • shatter head was chosen so that the pellet would penetrate into the plasma core even at larger displacements, as the head points into the direction of the movement.The timing of the SPI trigger provided reproducible results with constant delays between the triggering and the shattering.
Results.An overview of the current evolution and axis motion is given in figure 2, where the time of the first light measured by the AXUV diodes is marked with a dashed and the time of the TQ with a dotted line.The VDE motion is stopped by the TQ as shown in 2(b).Note that the current centroid is calculated from function parametrization [39], so that it can be influenced by large amounts of halo current that occur during the disruption.
The stop of the movement reduces the vertical moment change of the disruption and thus the force as shown in figure 3. # 40957 deviates from the linear regression, as the maximum of ∆M IZ is determined by the I p spike and not the displacement.As the edge safety factor increases due to the fast decay of the plasma current, the poloidal halo currents are greatly reduced as shown in figure 3(b).
In previous hot VDE experiments, the halo current was mainly observed in the divertor [38] as shown in the first row of figure 4.During the Neon mitigated experiments, however, significant halo current can also be measured in the heat shield (second row of figure 4).This is a strong indication that the halo current width is widened during mitigated disruptions as already observed in DIII-D [33].
The electron density n e rises up to a value of around 5 × 10 20 m −3 as measured by the CO2-interferometry [28,40].The radiated energy during the TQ and CQ is 0.6 MJ-0.7 MJ as measured by the bolometry system [41] makes up around 50% of the total energy [38].

Simulations and comparison to experiments
On the basis of the experiments in AUG and JET discharges described in [34], we performed MHD simulations for ASDEX Upgrade (AUG), JET and ITER.In the following, the MHD model is described, followed by the simulation setup and the results for JET and AUG.

Model and parameters
The MHD model.JOREK [42] is an extended MHD code that has been coupled to the EM code STARWALL [43,44] to take into account the effects of conducting structures around the plasma and to do free-boundary simulations [12].The reduced MHD version of JOREK, including effects of impurities [45], is used for the simulations presented here.Parallel flows have been neglected for numerical reasons.Also, two temperature effects are not included, so that in the following T, designates the total temperature (T = T e + T i with T i = T e ).The code has been used and benchmarked before for studies of VDEs in 2D [46] as well as 3D [37,38,47].Note that the convention of the sign of I p is machine dependent and can differ from the JOREK coordinate system.As the sign does not play a role for the analysis, we will only consider the absolute value in the following and likewise only the absolute value of I p ∆z.
Halo model and boundary conditions.Halo currents play a major role in the VDE dynamics and the resulting wall forces.In our model, the halo current is determined by the parallel electric field and the plasma resistivity in the halo region through Ohm's law.Thus, the temperature in the halo region plays a determining role.The halo temperature is given by the balance between Ohmic heating, impurity radiation and parallel losses along open field lines as described by Kiramov et al [48].Since we do not consider parallel flows in the work at hand and the impurity density in the halo region is experimentally difficult to assess, we study the influence of the halo temperature by modifying the parallel losses with two different setups.
1. Dirichlet boundary conditions and scan on the parallel conductivity (κ ∥ ).In this setup, the temperature at the wall is set to a low value (≈1 eV) and the parallel losses are determined by the parallel conductivity (κ ∥ ).The parallel conductivity model corresponds to the Spitzer-Haerm conductivity and the scans are performed by modifying the minimum (electron) temperature T min entering into the model, so that Therefore, a larger T min corresponds to a higher parallel heat conductivity in the SOL and thus, to larger conductive losses.2. Sheath boundary conditions and scan on the wall electron density.In this case, we impose a boundary condition on the heat flux with equation (7).To avoid the development of large gradients at the boundary, κ ∥ is set to a large constant value that leads to the homogenization of T along the field lines.Thus, the heat flux normal to the boundary q • n is determined by the field line temperature and the density at the wall: where γ sh is the sheath transmission coefficient [49], n wall e the electron density at the wall and c s the sound speed √ γk B T/m i .Note that we do not simulate self-consistently the evolution of the density at the wall, which would require to include neutral transport and recycling, but we use Dirichlet conditions instead.To study the influence of n wall e on the boundary condition, we scan γ sh (which is equivalent to a scan of n wall e at fixed γ sh ).
As found in the COMPASS tokamak [50], the halo current may also be constrained by the ion saturation current (or particle flux) at the plasma-wall interface if the ion density is not sufficiently large.Constraining the current density with the particle flux and thus, a fully self-consistent modeling of the halo current density would require sheath boundary conditions with the challenges described above which we leave for future work.

Simulation setup
The simulation setup for the different devices is described in the following.

Setup for ASDEX Upgrade simulations
Conductors and geometry.A polar grid with 90 radial and 120 polar elements is chosen to approximately fit the area inside the plasma facing components.The PF coil system, the VV and the passive stabilizing loop (PSL) are modeled [51].A thin-wall resistivity η w /d w of 30 µΩ has been used, which is smaller than the stainless-steel value of 50 µΩ to take into account the characteristics of the VV, resulting in a resistive wall time of several ms [52].Similar as in JET, the toroidal resistance is increased by bellows in between the sectors [53], so that the induced eddy current by the VDEs close inside the sector and thus have a larger effective path.It has been shown [38], that the movement of the hot VDE can be reproduced with the simulation and the current in the PSL evolves consistently with the experimental measurements.
Equilibrium and setup.The setup of the simulations and the equilibrium is the same as the 800 kA, 2 T discharge # 39655 described in [38].It is an Ohmic L-mode plasma with moderate elongation of κ = 1.8.To closely match experiments, the VDE is initiated at 2.97 s by a forced downward displacement of 2 cm in 20 ms by the virtual position controller in JOREK, while the other PF coils are frozen.Afterwards, also the current in the control coil is constant.

Setup for JET simulations
Conductors and geometrical setup.Similar to the AUG case, the plasma is modeled with a polar grid.The grid boundary does not fully correspond to the first wall contour in JET (see figure 5), but we consider it to be close enough for the main purposes of this paper.For the JET VV, we employ a simple axisymmetric thin wall model with an effective resistivity (η w ) and thickness (d w ) such that η w /d w = 430 µΩ [54], resulting in a resistive wall time of 5 ms [55].Since ferromagnetic conductors can presently not be modeled in JOREK-STARWALL, the magnetic effect of the iron core in JET is not taken into account, which could lead to potential differences with the experiments.Due to the absence of the iron core, the coil currents were adapted to optimally match the EFIT equilibrium [12].Initial equilibrium.All the JET runs shown in this paper are based on an EFIT reconstruction of the JET shot #95108 with I p = 1.1 MA and B = 1.2 T at t = 64.018s.Parabolic profiles for pressure and toroidal current density were chosen for the EFIT equilibrium reconstruction and thus detailed kinetic profiles are not taken into account.However, these will not play a role during the CQ phase since we introduce an artificial TQ and current profile flattening when impurities are injected.A spatially constant electron density of 1.5 × 10 19 m −3 is chosen throughout the domain to match the measured average density.The poloidally averaged profiles for temperature and toroidal current density as well as the safety factor are shown in figure 5(b).coil system and the passive components as described in [56], where these components were benchmarked with the DINA code.
Equilibrium and setup.The initial equilibrium for the presented simulations is an ITER 15 MA L-mode reference modelled with CORSICA [57] (H26 case) and used for SPI studies [58].All the ITER simulations are run with χ ⊥ = 1.3 m 2 s −1 .The upward VDE is triggered by a current perturbation added to the in-vessel coils.

Simulation phases
The simulations presented in this manuscript consist of the following phases Initial VDE.During this phase, the plasma drifts vertically while keeping its thermal and magnetic energy.The initial motion may be guided by modifying the PF coil currents or by slight profile changes.The parameters for this phase are the ones indicated in the previous section.

Artificial TQ.
In experiments, the injection of Neon leads to a quick thermal collapse.Such an effect is reproduced in the simulations by increasing the thermal diffusivities in order to obtain a desired TQ time.To facilitate the temperature drop, the Ohmic heating term is switched off during this phase.Once the electron temperature reaches values in the order of 10 eV in the plasma core, the thermal diffusivities are set back to their pre-TQ values and this phase concludes.
Current profile flattening and impurity ramp up.The experimental I p -spike is reproduced by flattening the current density profile with a helicity conserving term [59,60].The radial extension and strength of the flattening are free parameters that are tuned to match the experimental I p -spike.Once the TQ is complete, the impurity density is ramped up in the full simulation domain until the desired Neon density is obtained.
The ramp-up must be fast enough to ensure a balance between Ohmic heating and impurity radiation.Otherwise Ohmic reheating would lead to temperatures of hundreds of eVs, resulting in slow I p decay rates.The plasma temperature and thus the CQ decay rate is determined by the chosen impurity density.The distribution of the impurities is assumed to be uniform and not relevant as long as they can cool down the core sufficiently.In the absence of experimental measurements, the Neon density is chosen to match the initial CQ rate right after the I p -spike (e.g.see figure 7).Note that the impurity ramp-up and the current flattening are not simulated together with the TQ in order to facilitate the simulations.Otherwise, non-physical current distributions (e.g.thin current sheets) can occur in 2D simulations when the edge undergoes a radiative collapse while the core remains dominated by Ohmic heating.
CQ phase.After the chosen Neon density is attained, the impurity and main ion densities are kept constant and the artificial current profile flattening is switched off.The parameters during this phase are the same as in the initial VDE phase.However, scans of different parameters regarding the SOL properties and boundary conditions are performed during this phase as will be shown later.
For unmitigated VDEs, we do not apply most of these phases.Since the simulations are axisymmetric, the model has no restrictions on the q 95 evolution by MHD instabilities, and thus, unmitigated VDEs may reach values well below unity.In reality, 3D MHD activity will flatten the current profile, restricting the minimum q 95 [38].Only for the unmitigated VDE simulations, we apply a spatially constant hyperresistivity whenever q 95 drops below a given threshold (1 for JET) to mimic this effect.

The role of the halo temperature and width
As described in section 4.1, the role of the heat flux is studied by means of two different boundary conditions described in the following.Furthermore, the role of the halo current width is investigated.The simulations here follow the steps described in 4.3.

Minimum χ ∥ .
In the first case, Dirichlet boundary conditions are used for the density and temperature, while the minimum parallel heat conduction is varied.With larger parallel conduction losses (i.e.larger T min ), the halo temperature is reduced and consequently its resistance increases.Such an increase in resistance accelerates the halo current decay, which causes a faster global I p decay and a faster vertical motion since halo currents stabilize the VDE.
The latter can be observed in figure 8(a), where the influence of the parallel losses is shown for the JET case.When the parallel losses are not important, the current centroid stagnates as described in section 2. The edge safety factor then increases due to the fast CQ and the stagnation of the current centroid (i.e. the core current density does not increase due to the large halo currents).When the heat conduction becomes dominant (e.g. at T min = 100 eV), the plasma behaves similar to a hot VDE.In this case, the decreasing q 95 leads to larger poloidal halo currents.That explains the counter-intuitive observation that the poloidal currents are larger for the cases with smaller SOL temperature.
Figure 8(b) shows how the current moment, the vertical force and the poloidal halo currents behave with the varying parallel conduction.The cases with higher conduction (higher T min ) show a stronger movement, and thus a larger current moment.This increases the vertical force on the structures.Nevertheless, the current moment is significantly smaller than in unmitigated VDEs (see section 4.5) due to the current decay during the motion.
The same scan has been performed based on the AUG discharge # 40955, where T min was varied in the range 5 eV-35 eV. Figure 9(a) shows the influence of the parallel conduction on the R and Z position of the current centroid, the plasma current and q 95 .Note that the experimental values for q 95 and the current centroid position are only indicative as they result from fitting magnetic measurements to pre-calculated equilibria and thus are not fully reliable during disruptions.Thus, no quantitative comparisons can be done here.However, the figure indicates that the current centroid position first drops as a result of the TQ and then stagnates during the CQ phase, while the axis (shown with dashed lines in the same plot) continues to move downwards.The experimental CQ time is recovered with the lowest T min value, thus indicating that the energy lost by parallel conduction is negligible during the CQ.Due to the fast I p decay q 95 increases after the injection of impurities.Since the experiment only allows for an estimate of vertical forces [38], the comparison can only give an indication that the order of magnitude is comparable.Similar to what was already observed in figure 8  would lead to induction of more halo currents.On the other hand, the faster movement at lower T SOL reduces q 95 , which increases the poloidal halo currents.Figure 10 shows the radiation power in JOREK and AUG as well as the line integrated density.Due to the simulation setup, only the radiation during the CQ phase can be compared, where peak and duration can be well reproduced.As an Ohmic L-mode discharge, however, the magnetic energy poses the majority of the total stored energy, so that the CQ should contribute most to the radiated energy.However, there is a discrepancy in the line integrated electron density.

Boundary density scan.
In the second scan, the Dirichlet boundary conditions are replaced by sheath boundary conditions on the temperature as described in section 4.1.
Here, we study the effect of the density at the wall n wall e on the parallel losses.As it can be shown in the n wall e for JET in figure 11, similar conclusions are drawn as for the cases with Dirichlet conditions.Larger parallel losses lead to faster dynamics, less halo current stabilization and larger current moments and forces.Assuming a typical value of γ sh of 8, the simulations indicate the electron densities at the wall may be in the range of 10 20 -10 21 m −3 which are similar to the ones estimated for COMPASS unmitigated VDEs [61].
AUG results for the same scan are shown in figure 12.The results are comparable to the ones obtained with the scan of the parallel conduction.Again, the case with a low boundary heat flux corresponding to a density of 1.25 × 10 20 m −3 results in the experimental CQ time.Additionally, the evolution of the vertical moment is plotted, which shows that it is reduced with density, as the current centroid moves less in this case.Predictive results of ITER simulations in a scan against the value of the electron density at the wall (n wall e ) are shown in figure 13.Also for this scenario, we can see that a larger heat flux at the boundary leads to higher vertical forces, current centroid displacement and I p ∆Z curr .

Scan in the halo width.
As a consequence of the applied boundary conditions that do not allow for a selfconsistent density and temperature in the SOL, also the halo current width is not self-consistent.This has been observed before when comparing AUG simulations to the experiment [38], where the halo width was overestimated in the JOREK simulations.To show the influence of the halo current width in the simulation, a space dependent resistivity is introduced in the AUG simulations.For this, two limiting points are chosen beyond which the resistivity is increased by a factor of 100 to suppress halo currents far from the contact point with the structures.In practice, ψ is calculated at those points to identify the corresponding flux surface and the resistivity η is smoothly increased when crossing this surface.
Limiting the halo current width to a narrow region, corresponding only to the divertor or less, changes the dynamics significantly.The current centroid and the magnetic axis then move at a similar rate as the toroidal halo currents decay in the outer regions as shown in figure 14(a).The motion accelerates and q 95 decreases due to the shrinking core area, enhancing the poloidal halo currents.As the current centroid continues its trajectory, ∆M IZ increases further 14(c) and likewise the vertical force.
However, figure 15 shows that without restricting the halo width, the normal halo current to the boundary is comparable in its magnitude and extent to experimental shunts measurements in tiles on the heat shield, the divertor and on the PSL.While the halo current width in hot VDEs is relatively narrow as shown in section 3, the halo width is broader in the mitigated cases, reaching up to the heat shield.For a narrow halo region, the current density is larger than in the experiment and also the CQ time is reduced due to the faster movement and opening of the flux surfaces.

Global force reduction and influence of the injection delay
In this section, we study the influence of the delay in the pellet arrival on the global wall force.For that purpose, we inject the impurities at different time points (or vertical positions) during a reference unmitigated VDE.
JET.In figure 16(a), we show the simulated cases together with the available experimental data for the discharges described in by Gerasimov et al [34].The impurities are injected in the simulations according to the experimental times in which the shards reach the plasma.The setup with sheath boundary conditions, γ sh = 8 and n wall e = 7.5 × 10 19 m −3 was chosen for the 3 mitigated simulations corresponding to the discharges #95108, #95109 and #95153.Regardless of the  AUG results: when the halo width is restricted to a narrow region, the current centroid (a) continues its movement as the toroidal halo current decays leading to an increase of M IZ (c).As q 95 decreases the poloidal halo currents (b) are larger in this case.injection timing, the same behavior is observed for the mitigated cases.When the impurities arrive at the plasma, the current centroid stagnates, I p starts to decay, q 95 increases and |∆M Z | decreases.The maximum wall force is proportional to the maximum current moment (figure 16(b)) in agreement with theoretical considerations (section 2).Interestingly, the CQ time and the maximum poloidal halo current do not show a clear correlation with the maximum force.We conclude that for the JET case, F Z,max is determined by the vertical current moment at the time of the thermal collapse.In the case of the unmitigated VDE (discharge #19110), the maximum vertical moment |∆M Z | and its corresponding maximal force are determined by the minimum q 95 (since fixing q 95 implies limiting Z curr ).The similarity of M Z with experimental values indicates that the experimental lower threshold of q 95 is around unity as previously thought [36].ASDEX Upgrade. Figure 17(a) shows the time evolution of the SPI mitigated discharges (# 40955-40957) and the hot VDE (# 39655) as well as a case, where the TQ was initiated without a prior displacement (only simulation).These cases were carried out with sheath boundary conditions The evolution of plasma current, q 95 , the vertical current centroid position, the change of the current moment and the vertical force.The current centroid positions are qualitatively similar, but the excursion in the simulation is slightly larger.(b) The empty squares represent the scans of # 40955 in the SOL temperature for the current quench time, the maximum poloidal and toroidal halo current and the maximum of the vertical force against the maximum change int the current moment.The poloidal halo current magnitude is significantly reduced in the mitigated disruptions.As for JET, we find a good agreement with experimental data and a linear dependence between |Ip∆Zcurr| max and F Z,max .corresponding to n wall e of 1.4 × 10 20 m −3 .The current centroid evolution can be reproduced well before the injection of the SPI, while they only show a similar behavior during the CQ phase.Note, however, that the experimental values cannot be determined exactly.The poloidal halo currents are comparable for all mitigated discharges independent of the injection time, while the values for the hot VDE are larger because q 95 has a low value ≈2 for this case [38].Due to the uncertainty of Z curr in the experiment, also the current moment cannot be determined exactly and the curves during the CQ phase deviate from the simulation results.
Figure 17(b) shows the simulation results of the maximum vertical force and halo current for the cases with different injection times as well as the hot VDE.The squares indicate the scan of the parallel conduction and boundary density for # 40955.The force is proportional to the maximum of the current moment.On the other hand, the poloidal halo currents have comparable magnitudes for all the cases, meaning that the vertical force cannot mainly depend on them.

Predictions for ITER.
In figure 18, we show the results of simulations of the ITER scenario described in section 4.2.3.
Here, sheath boundary conditions have been considered, with γ sh = 8 and n wall e = 6 × 10 20 m −3 .From the top to the bottom of figure 18(a), we show the time traces of the plasma current, the edge safety factor, the vertical displacement of the current centroid, I p ∆Z curr and the vertical wall force.The solid dark blue line represents the dynamics of an unmitigated VDE.The other colored lines correspond to the time traces of simulations where the impurities have injected at approximately 90, 120, 230 and 300 ms respectively, to different displacements of current centroid and magnetic axis.As for the AUG and JET simulations, these mitigated cases present the same common features.Once the impurities reach the plasma, the plasma current decays, q 95 increases and I p ∆Z curr decreases.The vertical displacement of the current centroid modifies the slope of its trajectory until it reaches a maximum and then it begins to decrease.Similar to the JET simulations, the current diffusion was tuned by the hyperresistivity to prevent q 95 from decreasing below unity in the unmitigated case.
In figure 18(b), the maximum of the vertical force measured in JOREK simulations (colored triangles) is shown as a function of the maximum value of I p ∆Z curr .Note that the triangles have the same colors of the corresponding time traces displaced in figure 18(a).Unlike the results shown for JET and AUG, the maximum of the vertical force is not proportional to the maximum of I p ∆Z curr due to the shielding of the wall forces described by the second term in equation (1).In JET and AUG, the decay time of wall currents is smaller than the time of the disruption τ w ≪ τ D = τ VDE + τ CQ .Therefore, the shielding is not important.However, in the ITER cases, the disruption time transitions from τ D ≪ τ w with τ w ≈ 235 ms [31] for an early mitigation to τ D ≈ τ w in the hot VDE case.For an early injection, the shielding from the wall is more important than for a late injection or the hot VDE.To show that this is indeed the result of the shielding, we apply the formula given by Miyamoto [31, equation (6)] that estimates the force reduction of F p,c by eddy currents in the wall.This formula was applied to I p ∆Z curr in order to obtain the factor of proportionality α to F p,c .Figure 18(b) shows the theoretical prediction of the vertical force, where equation ( 6) from [31] was applied to αI p ∆Z curr and compared to the values obtained from the simulation.From this, it becomes evident that the shielding is responsible for the force reduction, which is more than linear for early injections in ITER due to the shielding effects of a conductive wall.For an early injection, the maximum wall force is smaller than 5 MN.The effect of the conductive wall on the force reduction has also been described here [31,62], where a decrease of the vertical force with τ CQ was found.The reduction of the wall force in the limit of τ w ≪ τ D was also predicted by Pustovitov et al [63].

Conclusion
In this manuscript, we proposed a new mechanism explaining the force reduction during mitigated disruptions based on the reduction of the vertical current moment.Simulations with the MHD code JOREK confirm that the vertical moment is reduced due to large parallel currents in the halo region that stop the movement of the current centroid as previously described in [33].While the plasma core continues its vertical motion in accordance with theoretical predictions, the geometric center of the currents remains stationary after the injection of impurities.This is supported by experiments in AUG and JET, where a stationary current centroid is observed after an SPI.Even if the plasma is already significantly displaced during a VDE, the motion stops shortly after the injection.As the vertical force is proportional to I p ∆Z curr , this explains the force reduction in the experiment.Another confirmation of the theory is the widening of the halo width in mitigated discharges in comparison to hot VDEs that is observed in AUG and DIII-D [33], which was not investigated in JET.
While the simulations were only carried out in 2D and with a uniform impurity density, the free parameters (impurity density and heat flux) were adapted to include the necessary features to explain the experimental behavior.The impurity density was chosen to fit the initial CQ rate in JET and AUG.The choice is confirmed by measurements of the radiated power during the CQ in AUG by foil bolometers, which are in the same order as in JOREK.The latter parameter was scanned by varying either the parallel conduction in the SOL or the heat flux on the boundary.In disruptions mitigated by a large quantity of impurities, the main energy sink is supposed to be radiation [23].Therefore, a good match with the experiment is obtained for low values of the conducted energy to the boundary.Altogether, the current centroid motion, the CQ time and the magnitude of the vertical force is qualitatively in agreement with the experiment.In JET, the measurement of the vertical moment can be directly compared to the JOREK simulations for different injection times with a good correspondence.In both AUG and JET, it was confirmed with JOREK that the vertical force scales linearly with the maximum of the current moment and thus are reduced the earlier the injection is done during a VDE.
Predictive simulations for ITER have shown that the vertical force can be limited to 5 MN, if an early enough SPI occurs.For the ITER cases, the shielding by the induced eddy currents in the wall can reduce the total vertical force since the decay time of the wall currents is longer than the disruption time.
While these results provide a hopeful outlook on future experiments as the mitigation of the vertical force is now better understood, there are still open questions to answer.First, the reduction of the vertical force relies on an effective injection of impurities to reduce the core temperature to the order of tens of eV, which prevents the re-induction of current in the core.We do not model the impurity injection and assimilation realistically to simplify the simulation.However, present experiments show a reliable control of the CQ time with MGI and SPI.Second, while the poloidal halo currents are significantly reduced due to an increased q 95 , it has to be verified that the heat loads and the local J × B forces are sufficiently mitigated.Finally, this study leaves out RE currents.In the case of a focused RE beam of several MA moving vertically, ∆M IZ is not reduced and thus the vertical force could still be significant.Estimates for ITER show a maximum RE current of 10 MA [64] which would lower the force to 2/3 of the unmitigated value.However, this is likely to be an upper estimate.RE termination events at low q 95 could lead to the reformation of an Ohmic current and the transition to a thermal resistive plasma where forces are well mitigated.Also, the presence of toroidal halo currents during the CQ could be beneficial for RE mitigation since they could reduce the available magnetic energy for avalanching in the confined region.Studies to investigate the impact of an RE beam are left for future work.

Figure 1 .
Figure 1.The TQ is initiated at 44 ms and the temperature profile is flattened into the SOL, so that a large part of the current is carried by halo currents.(a) The magnetic axis (blue, dashed) moves monotonically downwards, while the current centroid (red) stays closer to the midplane and reverses after the LCFS is lost.The toroidal halo current I h,ϕ can make up a large part of the initial plasma current.(b) The toroidal current distribution of an AUG simulation during the CQ phase.The trajectory of the magnetic axis and current centroid are shown in white and red respectively.

Figure 2 .
Figure 2. Based on the hot VDE AUG discharge # 39655, a TQ was triggered intentionally by a 100% Neon SPI at three different times during the displacement.The time of the pellet arrival (first AXUV light) is shown by a dashed line and the dotted line indicates the time of the TQ onset.(a) Ip decays soon after the pellet arrives and the halo currents are reduced with early injections.(b) The current centroid does not move significantly after the injection in comparison with the hot VDE # 39655.

Figure 3 .
Figure 3. (a) The forces increase with the change of the vertical moment Ip∆Z, where a linear fit is indicated.(b) The sooner the SPI is injected the more the halo fraction (HF) is reduced.

Figure 4 .
Figure 4.The halo current as measured by the shunt diagnostics in AUG is shown here, where the time is shifted with respect to the start of the VDE at t 0 .The normal current to the PFCs is wider in the SPI mitigated discharges (# 40955-# 40957) compared to previous hot VDE experiments (# 39655-# 39722).Note, that some parts of the heat shield are not covered with diagnostics and the value was interpolated for the SPI mitigated cases.

Figure 5 .
Figure 5. (a) Geometrical setup for the JET simulations.The JOREK grid is colored in green, the red line denotes first wall contour and the poloidal contour of the employed axisymmetric thin vacuum vessel model is depicted yellow.The separatrix of the initial equilibrium is shown for reference.(b) Flux surface averaged equilibrium profiles of temperature, current density and safety factor as a function of the normalized poloidal flux.

Figure 6 .
Figure 6.Geometrical setup for the ITER simulations.The JOREK flux surface aligned grid is shown in green, the ITER first wall in red, while the yellow line represents the vacuum vessel (VV).The black line indicates the separatrix of the initial equilibrium.

4. 2 . 3 .
Setup for ITER simulationsConductors and geometrical setup.Unlike in the previous two cases, a flux surface aligned grid presented in green in figure6.It has 67 radial and 128 poloidal grid elements inside the separatrix.The VV has been modelled together with the PF

Figure 7 .
Figure 7. Reproduction of experimental Ip-spike and initial CQ rate in the JET simulations.The solid lines correspond to simulations with different T min (in eVs) used for the Spitzer-Haerm parallel conduction κ ∥ (max(T, T min )) and the dashed line to the experiment.The choice of a colder halo (larger T min ) does not influence the initial Ip decay rate, which is given by the chosen of Neon density.

Figure 8 .
Figure 8. JET results: scan in the parallel conduction and comparison to shot #95108.(a) From top to bottom: time traces of the Z and R displacements of the current centroid, plasma current and edge safety factor.The current centroid does not move significantly after the TQ, and the values reproduce the main experimental features like the CQ time and ∆Zcurr.(b) From top to bottom: time traces of the toroidal and poloidal halo current, change of the vertical current moment and global vertical force on the vacuum vessel.The solid lines correspond to simulations with different T min used for the Spitzer-Haerm parallel conduction κ ∥ (max(T, T min )) and the dashed line to the experiment.The vertical force varies with M Z , but does not depend strongly on the poloidal halo current.
(b), figure 9(b) shows how the toroidal halo currents do not vary much with different injectimes prior to the loss of the LCFS, as a faster movement

Figure 9 .
Figure 9. AUG results: scan in the parallel conduction and comparison to shot #40955.(a) From top to bottom: current centroid Z position and R displacement, plasma current and edge safety factor.The current centroid does not move significantly after the TQ, and the values reproduce the main experimental features like the CQ time and ∆Zcurr.(b) From top to bottom: time traces of the toroidal and poloidal halo current, change of the vertical current moment and global vertical force on the vacuum vessel.

Figure 10 .
Figure 10.Evolution of the radiation obtained from bolometer signals and the line integrated density obtained from the vertical interferometer lines.The Neon impurity density in case #40955 was chosen to fit the initial CQ time.Also, the peak of the experimentally measured radiation is comparable, while the line integrated electron density in the simulation is lower.

Figure 11 .
Figure 11.JET results: scan on the electron density at the wall to vary the boundary heat flux and comparison to experiments.(a) From top to bottom: time traces of the Z and R displacements of the current centroid, plasma current and edge safety factor.(b) From top to bottom: time traces of the toroidal and poloidal halo current, change of the vertical current moment and global vertical force within the vacuum vessel.The solid lines correspond to simulations with different electron densities used in the heat flux boundary condition (with γ sh = 8) and the dashed line to the experiment.The results are comparable to the scan of T min , κ ∥ , therefore the choice of the BC do not have a large influence.

Figure 12 .
Figure 12.AUG results: scan in the n wall e and comparison to the discharge #40955.(a) From top to bottom: time traces of the Z position and R displacements of the current centroid, plasma current and edge safety factor.(b) From top to bottom: time traces of the toroidal and poloidal halo current, change of the vertical current moment and global vertical force within the vacuum vessel.The results here are comparable with the scans in the minimum parallel heat diffusion, showing a slower movement with larger SOL temperature.Again, we can observe that the current centroid movement is limited after the TQ.

Figure 13 .
Figure 13.ITER results: scan in the electron density at the wall (n wall e ) in simulations of a VDE ITER, where the TQ was triggered at a displacement of 40 cm.(a) From top to bottom: time traces of the vertical and horizontal displacements of the current centroid, total plasma current and edge safety factor.(b) From top to bottom: time traces of the toroidal and poloidal halo currents, time traces of the current moment variation and the vertical force.The results agree with the AUG and JET simulations, where the current centroid remains stationary after the injection and the force is reduced.

Figure 14 .
Figure14.AUG results: when the halo width is restricted to a narrow region, the current centroid (a) continues its movement as the toroidal halo current decays leading to an increase of M IZ (c).As q 95 decreases the poloidal halo currents (b) are larger in this case.

Figure 15 .
Figure 15.AUG results: the experimental halo width of the AUG discharge # 40957 is compared to the normal halo current to the JOREK computational boundary.The edges of the divertor is indicated by dashed lines.The shunt measurements on the heat shield (on the HFS) are interpolated for the regions where no measurements are present.Without restricting the width, the halo current magnitude and extent are comparable.

Figure 16 .= 7 . 5 ×
Figure 16.JET results: reduction of the global vertical force in JET by injecting Neon at different vertical displacements.(a) From top to bottom: time traces of the plasma current, edge safety factor, vertical displacement of the current centroid, the change of the vertical current moment and the total vertical wall force.Solid lines indicate simulations, while dashed lines indicate experimental results.The simulations were performed with the sheath boundary condition approach with n wall e = 7.5 × 10 19 m −3 .A good agreement with the experiment for various injection times is found.(b) From top to bottom: CQ time (20-80 definition), the maximum toroidal and poloidal halo current and the maximum vertical force as a function on the maximum change in vertical moment.The shaded area corresponds to the scans performed for case #95108 in figures 8 and 11.The three black triangles indicate the simulations of shots #95109, #95135 and #95110.It can be seen that F Z,max is proportional to |Ip∆Zcurr| max , while it is independent of the amount of poloidal halo currents or τ CQ here.

Figure 17 .
Figure17.AUG results: comparison to the experimental results for multiple AUG shots.(a) The evolution of plasma current, q 95 , the vertical current centroid position, the change of the current moment and the vertical force.The current centroid positions are qualitatively similar, but the excursion in the simulation is slightly larger.(b) The empty squares represent the scans of # 40955 in the SOL temperature for the current quench time, the maximum poloidal and toroidal halo current and the maximum of the vertical force against the maximum change int the current moment.The poloidal halo current magnitude is significantly reduced in the mitigated disruptions.As for JET, we find a good agreement with experimental data and a linear dependence between |Ip∆Zcurr| max and F Z,max .

Figure 18 .
Figure 18.ITER results: (a) from top to bottom: time traces of the total plasma current, edge safety factor, vertical displacement of the current centroid, Ip ∆Zcurr and total vertical wall force.The time traces of an unmitigated VDE (UVDE), are compared with those of simulations where impurities have been injected at different times.(b) Maximum vertical force as a function of the maximum change in Ip ∆Zcurr.Note that the triangles have the same colors of the corresponding time traces shown in the left plots.These results, obtained from JOREK simulations, are compared with the theoretical predictions taking into account the vessel shielding effects.