1. Introduction
One of the main issues threatening the success of future reactor-scale tokamaks is disruptions. It is the sudden loss of confinement, where the plasma rapidly ejects a large portion of its energy content onto the first wall, which exposes the device to excessive mechanical stresses, heat loads and can lead to the formation of a runaway electron beam (Hender et al. Reference Hender2007; Breizman et al. Reference Breizman, Aleynikov, Hollmann and Lehnen2019). Unmitigated disruptions could potentially cause severe damage to the device, and modelling such events is thus crucial to assess the effectiveness of various mitigation techniques and to optimise termination scenarios.
Shattered pellet injection (SPI) employs hydrogen isotopes (protium or deuterium) with minority neon admixture, frozen into cryogenic pellets, which are launched into the vessel at speeds of several hundred metres per second. Upon entry, the singular pellet is broken into smaller pieces by a shattering unit to increase pellet material ablation and assimilation into the plasma. The exact pellet composition, penetration velocity and the size distribution of pellet fragments has a major impact on the material assimilation (Jachmich et al. Reference Jachmich2023) and disruption mitigation efficiency (Heinrich et al. Reference Heinrich2025; Patel et al. Reference Patel2025; Tang et al. Reference Tang2025). Modelling is crucial to understand the complex physics and to optimise the termination scenarios.
The ITER baseline disruption mitigation system will be based on SPI due to its ability to promptly inject material into the plasma to thermally radiate its stored energy, reduce electromagnetic loads on the surrounding conducting structures and to inhibit the generation of runaway electrons (Hollmann et al. Reference Hollmann2014; Jachmich et al. Reference Jachmich2021; Lehnen et al. Reference Lehnen2015; Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024a ). As the mitigation of the various disruption consequences, however, poses conflicting criteria, finding an optimised scheme is not straightforward, requiring experimental and theoretical input. One such conflict could be if different SPI penetration speeds and fragment distributions were required for mitigation thermal loads, vessel forces or runaway electrons. In support of its development, a highly flexible SPI system was installed at the ASDEX Upgrade (AUG) tokamak for the 2021/2022 SPI experimental campaign to investigate the impact of different fragment size and velocity distributions on the disruption dynamics (see Dibon et al. Reference Dibon2023; Heinrich et al. Reference Heinrich, Papp, de Marné, Dibon, Jachmich, Lehnen, Peherstorfer and Vinyar2024). The AUG system has so far been using deuterium as the primary hydrogen isotope. We do not expect major differences in the disruption dynamics between protium and deuterium SPI from the atomic physics point of view.
High-fidelity magnetohydrodynamic (MHD) simulations are routinely applied to disruption modelling, e.g. with Jorek (Nardon et al. Reference Nardon, Hu, Hoelzl and Bonfiglio2020; Hoelzl et al. Reference Hoelzl2021; Hu et al. Reference Hu, Nardon, Hoelzl, Wieschollek, Lehnen, Huijsmans, van Vugt and Kim2021; Tang et al. Reference Tang2025) or Nimrod (Kim et al. Reference Kim, Liu, Parks, Lao, Lehnen and Loarte2019; Sovinec et al. Reference Sovinec, Glasser, Gianakon, Barnes, Nebel, Kruger, Schnack, Plimpton, Tarditi and Chu2004), representing the high physics fidelity – high numerical cost – end of the spectrum. To carry out large parameter scans, aid design decisions and help focus expensive simulations where most needed, reduced models are also employed. In this paper, the one-dimensional fluid model of the Dream code (disruption runaway electron analysis model, by Hoppe, Embreus & Fülöp (Reference Hoppe, Embreus and Fülöp2021)) is used.
This code has been used for simulations of ITER (Pusztai, Hoppe & Vallhagen Reference Pusztai, Hoppe and Vallhagen2022; Lier et al. Reference Lier, Papp, Lauber, Pusztai, Särkimäki and Embreus2023; Ekmark et al. Reference Ekmark, Hoppe, Fülöp, Jansson, Antonsson, Vallhagen and Pusztai2024; Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024a
) as well as TCV (Wijkamp et al. Reference Wijkamp, Hoppe, Decker, Duval, Perek, Sheikh, Classen and Jaspers2023), JET (Järvinen et al. Reference Järvinen, Fülöp, Hirvijoki, Hoppe, Kit and Åström2022), DIII-D (Izzo et al. Reference Izzo, Pusztai, Särkimäki, Sundström, Garnier, Weisberg, Tinguely, Paz-Soldan, Granetz and Sweeney2022; Marini et al. Reference Marini, Hollmann, Tang, Herfindal, Shiraki, Wilcox, del Castillo-Negrete, Yang, Eidietis and Hoppe2024), STEP (Berger et al. Reference Berger, Pusztai, Newton, Hoppe, Vallhagen, Fil and Fülöp2022; Fil et al. Reference Fil, Henden, Newton, Hoppe and Vallhagen2024) and SPARC (Izzo et al. Reference Izzo, Pusztai, Särkimäki, Sundström, Garnier, Weisberg, Tinguely, Paz-Soldan, Granetz and Sweeney2022; Tinguely et al. Reference Tinguely2023). In this paper we demonstrate that Dream can explain experimental trends of major disruption characteristics at AUG, such as the evolution of the plasma current
$I_{\textrm {p}}$
or the radiated energy fraction
$f_{\textrm {rad}}$
.
The plume of fragments used as input for the modelling is generated by sampling fragment size and velocity distributions that depend on various injection parameters, such as the pellet neon fraction
$f_{\textrm {Ne}}$
and the injection speed
$v_{\textrm {inj}}$
, using the Parks model (Parks Reference Parks2016) as implemented by Gebhart et al. Reference Gebhart, Baylor and Meitner(2020a). Since pellet fragmentation is a statistical process, we simulate an ensemble of several realisations for each parameter set, to assess the impact that statistical variation in the fragment distribution has on the resulting disruption dynamics (sensitivity analysis).
The physics model used is detailed in §§ 2.1 and 2.2, followed by a description on the fragment sampling in § 2.3. An overview of the simulation set-up and the AUG SPI discharges selected for the experimental comparison is presented in § 2.4 and 2.5, respectively. A demonstrative scan in the pellet neon fraction is detailed in § 3.1, and the experimental comparison in § 3.2. The results and limitations of the model are discussed and concluded in § 4.
2. Methods
2.1. Physics model
A one-dimensional fluid model of the plasma available in the Dream code is applied. This is a finite-volume numerical simulation framework that has been developed for tokamak disruption modelling, with a particular focus on the runaway electron dynamics. It evolves the plasma current self-consistently, along with the thermal bulk of electrons, and charge states of the relevant ionic species – all in a realistic toroidally symmetric magnetic geometry. In this section, the relevant equations and physical assumptions are described in detail.
2.1.1. Plasma current
The plasma current density
$j_\parallel$
consists of an Ohmic component
$j_\Omega$
, and a component
$j_{\textrm {re}}$
carried by runaway electrons. Other non-inductively driven currents, such as the bootstrap current, are not included in this model as these are expected to convert to Ohmic current during a disruption. The total current density is related to the poloidal magnetic flux
$\psi _{\textrm {p}}$
, as described by Ampère’s law

where
$\mu_0$
is the permeability of free space,
$\langle \cdot \rangle$
denotes flux surface average,
$\boldsymbol{B}$
the magnetic field,
$R$
the major radius and
$V'$
is the spatial Jacobian (
$\textrm {d}{V}=V'\textrm {d}{r}$
being the volume between two infinitesimally close flux surfaces). These depend on the particular magnetic equilibrium
$\boldsymbol{B}(r, \theta )$
used in the simulation set-up, which is treated as static in this model. Note that radial coordinate
$r$
is defined as the distance on the midplane from the magnetic axis
$R={R_0}$
to the corresponding flux surface on the low-field side.
The poloidal magnetic flux
$\psi _{\textrm {p}}$
is evolved using a modified version of Faraday’s law of induction

This is the loop voltage, of which the first term on the right-hand side is proportional to the component of the electric field
$\boldsymbol{E}$
parallel to the magnetic field driving the Ohmic component of the current according to Ohm’s law
$j_\Omega /B=\sigma \langle \boldsymbol{E}\cdot \boldsymbol{B}\rangle /\langle B^2\rangle$
. The parallel electric conductivity is modelled as outlined in Redl et al. (Reference Redl, Angioni, Belli and Sauter2021). It is a neoclassical correction of the Spitzer conductivity
$\sigma \propto T_e^{3/2}$
accounting for trapped particle effects at arbitrary collisionality.
The second term in (2.2), derived by Boozer (Reference Boozer1986), involves a second-order derivative in the toroidal magnetic flux
$\psi _{\textrm {t}}=(2\pi )^{-1}\int V'\textrm {d}{r}\langle \boldsymbol{B}\cdot \boldsymbol{\nabla }{\varphi }\rangle$
, which is used as a radial coordinate. Given the relation in (2.1), this term is in fact a fourth-order radial derivative, and thus describes the hyperdiffusion of
$\psi _{\textrm {p}}$
. Whenever
$\Lambda$
is greater than zero, this term acts to relax
$j_\parallel /B$
while conserving the helicity content of the plasma. This leads to a transient increase in
$I_{\textrm {p}}$
as the current density radially redistributes in what is referred to as the ‘
$I_{\textrm {p}}$
-spike’. This emulates the effects that the break-up of magnetic field lines has on the current density, and is used to model the enhanced transport event during the thermal quench phase of a disruption with strong field stochastisation. This is explained in more detail in § 2.2.
As boundary condition for the plasma edge flux
$\psi _{\textrm {p}}(a)$
, we couple it to an ideally conducting wall at
$r=b\gt a$
. By doubly integrating (2.1) over the vacuum region
$a\lt r\lt b$
, in which
$j_\parallel =0$
, we can relate the plasma edge flux to that of the wall
$\psi (b)=\psi (a)+M_{\textrm {ab}}I_{\textrm {p}}$
, where
$M_{\textrm {ab}}=\mu _0{R_0}\ln (b/a)$
is the edge–wall mutual flux inductance taken in the cylindrical limit
$a/R\to 0$
. We consider a wall that is perfectly conducting, meaning no electric field is induced in the wall and thus the wall flux remains constant in time.
In addition to the Ohmic component, the plasma current can also contain a component driven by runaway electrons. Without resolving the electron distribution in momentum space, the runaway electrons are assumed to travel at the speed of light
$c$
, and their contribution to the total plasma current density is given by
$j_{\textrm {re}}=ecn_{\textrm {re}}$
, with
$e$
being the elementary charge. The evolution of the runaway density
$n_{\textrm {re}}$
also includes diffusive transport with the diffusion coefficient set to
$D_{\textrm {re}}=1\,\textrm{m}^2\,\textrm{s}^{-1}$
, and a set of source terms accounting for the runaway electron generation mechanisms. These include Dreicer generation, which is calculated by a neural network that has been trained on kinetic simulations (Hesslow et al. Reference Hesslow, Unnerfelt, Vallhagen, Embreus, Hoppe, Papp and Fülöp2019b
; Ekmark et al. Reference Ekmark, Hoppe, Fülöp, Jansson, Antonsson, Vallhagen and Pusztai2024); hot-tail generation, as calculated in § 4.2 of Svenningsson (Reference Svenningsson2020); and a model for the avalanche growth rate, which accounts for the effects of partially screened nuclei in non-ideal plasmas (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a
). A negligible generation of runaways is expected for the plasma scenarios considered in this work, but the respective terms are included for completeness as these models are readily available in Dream and are relatively inexpensive computationally. Runaway electrons were so far not observed in ASDEX Upgrade SPI experiments with neon and deuterium–neon injections. Our simulations help explain that this is due to insufficient generation, rather than due to runaway seed losses preventing avalanche.
2.1.2. Ions and SPI dynamics
Ions and neutral atoms are characterised both by a temperature
$T_i$
, and a set of densities
$n_i^{(j)}$
for every species
$i$
with atomic number
$Z_i$
, and charge state
$j=0,\ldots, Z_i$
. The densities are evolved according to

The first four terms describe ionisation and recombination processes, for which the coefficient
$\mathcal{I}_i^{(j)}$
denotes the rate of ionisation between the two charge states
$j\to j+1$
, and
$\mathcal{R}_i^{(j)}$
the rate of recombination between
$j\to j-1$
. These coefficients are extracted from the OpenADAS database (Summers Reference Summers2004).
The fifth term describes a source into the neutral state
$n_i^{(0)}$
of the injected material, as the fragments travel along their trajectories
$r_k(t),\, k=1,\ldots N_{\textrm {frag}}$
, depositing ablated material. The rate at which a fragment decreases in diameter
$d_k(t)$
is
$\dot {d}_k=-2\mathcal{G}_k\mathcal{M}/(\pi d_k^2 \rho N_{\textrm {A}})$
, where
$\mathcal{M}$
is the pellet molar mass,
$\rho$
the pellet mass density and
$N_{\textrm {A}}$
Avogadro’s number. The rate of the deposition of neutrals per unit time
$\mathcal{G}_k$
is calculated using a version of the neutral gas shielding (NGS) model, first derived by Parks and Turnbull (Reference Parks and Turnbull1978) and further extended upon by Parks (Reference Parks2017) to allow for mixed deuterium–neon pellets. The NGS model provides the scaling law
$\mathcal{G}_k\propto \lambda n_e^{1/3}T_e^{5/3}d_k^{4/3}$
, where
$\lambda$
(along with
$\mathcal{M}$
and
$\rho$
) is a function of the neon fraction of the pellet
$f_{\textrm {Ne}}=N_{\textrm {Ne}}/(N_{\mathrm{D_2}}+N_{\textrm {Ne}})$
.
The sixth and final term in (2.3) describes both advective and diffusive particle transport, for which the transport coefficients
$A_i$
and
$D_i$
are active during the thermal quench (similarly to
$\Lambda$
). More on this in § 2.2.
For each ion species, a temperature
$T_i$
, which is assumed to be same for all charge states, is evolved via collisional heat exchange with the bulk electron population and other ion species present in the plasma, according to

with the total density
${n_i}=\sum _{j=0}^{Z_i}n_i^{(j)}$
. Here,
$Q_{\alpha \beta }$
is the collisional heat exchange between Maxwellian species
$\alpha$
and
$\beta$
, which is is proportional to their temperature difference, as further detailed in Hoppe et al. (Reference Hoppe, Embreus and Fülöp2021).
2.1.3. Thermal electrons
The evolution of the electron temperature
$T_e$
is modelled using an energy balance equation

in which the first term on the right-hand side describes Ohmic heating, and the second radiative cooling due to line radiation, recombination radiation and bremsstrahlung. Integrating this term over the plasma volume yields the total radiated power, which we will denote as
$P_{\textrm {rad}}$
. Typically, this energy balance is dominated by Ohmic heating and radiative cooling, both of which are nonlinear functions of the electron temperature. The rates of energy loss
$L_i^{(j)}(n_e, T_e)$
for neon are retrieved from OpenADAS (Summers Reference Summers2004). Deuterium is instead assumed to be opaque to Lyman radiation, which reduces the radiative losses at temperatures below
$10\,\textrm {eV}$
; for this the loss rates are extracted from the AMJUEL database (Reiter Reference Reiter2020).
The third term,
$\mathcal{P}_{\textrm {ion}}$
, represents the change in potential energy due to the electron–ion binding energy removed by the electron bulk during ionisation and added during recombination. With the ionisation threshold energies
$\Delta W_i^{(j)}$
obtained from the NIST database (Kramida et al. Reference Kramida, Ralchenko and Reader2020), this contribution to the energy balance is calculated by adding up the terms
$\Delta W_i^{(j)}(\mathcal{R}_i^{(j)}-\mathcal{I}_i^{(j)})n_i^{(j)}n_e$
over all charge states.
Collisional heat transfer from the runaway electrons is accounted for with the fourth term, for which
$E_{\textrm {c}}=n_e e^3\ln {\Lambda }/(4\pi \varepsilon _0^2 m_e c^2)$
, where
$\ln\Lambda$
is the Coulomb logarithm,
$\varepsilon_0$
the permittivity of free space, and
$m_e$
the mass of the electron. The fifth term is the collisional heat transfer between electrons and ions
$Q_{ei}=-Q_{ie}$
, and the final term describes diffusive heat transport, with
$\chi _e$
being the heat diffusivity. A baseline value of
${\chi _e}=1\,\mathrm{m}^2\,\mathrm{s}^{-1}$
is used throughout.
2.2. Enhanced transport event
To emulate the effect of the magnetic field stochastisation during the thermal quench that causes the thermal energy loss, an enhanced radial transport of heat, particles and poloidal magnetic flux is applied at a time
$t_{\textrm {relax}}$
. As previously done in ITER studies by Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai(2024a), we take
$t_{\textrm {relax}}$
to be the time when the electron temperature drops below
$10\,\textrm {eV}$
somewhere inside of the resonant
$q=2$
surface. The choice of this particular criterion for triggering the enhanced transport event is motivated by the fact that the conductivity at
$10\,\textrm {eV}$
is sufficiently low to likely trigger a resistive MHD instability at the resonant
$q=2$
flux surface on a millisecond time scale.
In the simulation the transport event is emulated by increasing any transport coefficient
$\mathcal{Y} \in \{ \Lambda , A_i, D_i, {\chi _e} \}$
from their baseline value
$\mathcal{Y}_0$
up to some maximum value
$\mathcal{Y}_{\textrm {max}}$
, which then exponentially decays over a time scale of
$\tau _{\textrm {decay}}=1\,\textrm {ms}$
, according to

where
$\varTheta$
denotes the Heaviside function. In particular, this concerns the hyperdiffusivity
$\Lambda$
in (2.2), the ion and neutral particle advection and diffusion coefficients
$A_i$
and
$D_i$
in (2.3) and the electron heat diffusivity
$\chi _e$
in the energy balance (2.5). We omit the enhanced transport of runaway electrons, as all simulated cases considered here yield already negligible runaway currents. We have conducted a sensitivity study for the impact of
$\Lambda _{\textrm {max}}$
on current quench evolution in the range of
$\Lambda _{\textrm {max}}=[1\times 10^{-7}-1\times 10^{-3}]\,\mathrm{Wb^2}\, \textrm{m}^2\,\textrm{s}^{-1}$
. By selecting
$\Lambda _{\textrm {max}}=6.3\times 10^{-5}\,\mathrm{Wb^2}\,\textrm{m}^2\,\textrm{s}^{-1}$
, we observe an increase in the plasma current during the
$I_{\textrm {p}}$
-spike similar to the experiment (
$\Delta I_{\textrm {p}}\sim 20\,\textrm {kA}$
for the cases studied here). The impact of
$\Lambda _{\textrm {max}}$
in this range on the radiated energy fraction is between
$\pm 2\,\%$
and
$\pm 5\,\%$
, depending on the pellet composition. The maximum value of the other transport coefficients are set to
$A_{i,\textrm {max}}=-100\,\textrm{m}\,\textrm{s}^{-1}$
,
$D_{i,\textrm {max}}=100\,\textrm{m}^2\,\textrm{s}^{-1}$
and
$\chi _{e,\textrm {max}}=100\,\textrm{m}^2\,\textrm{s}^{-1}$
. Similar values were used in ASTRA simulations of massive gas injection scenarios in AUG by Linder et al. (Reference Linder, Fable, Jenko, Papp and Pautasso2020) to reproduce the evolution of line-integrated electron densities from experiment.
2.3. Fragment sampling
The generation of the fragment plume is performed in two steps: a rejection sampling of fragment sizes, utilising an analytical distribution function
$f(d)$
for the characteristic fragment diameter
$d$
following the pellet break-up; and then assigning velocities
$\boldsymbol{v}$
to the individual fragments. This results in a set of tuples
$\{(d_k, \boldsymbol{v}_k)\}$
, with
$k=1,\ldots, N_{\textrm {frag}}$
.
The fragment size probability distribution stems from a statistics-based model for the shattering of brittle materials, initially derived in the context of exploding munitions by Mott and Linfoot (Reference Mott and Linfoot1943), and then more recently by Parks (Reference Parks2016) in the context of shattered cryogenic pellets. This has subsequently been correlated with relevant injection parameters to experimentally measured fragment distributions to yield a modified analytical distribution for the fragment size (Gebhart et al. Reference Gebhart, Baylor and Meitner2020a
,Reference Gebhart, Baylor and Meitner
b
), which is used in this paper. For a cylindrical pellet of diameter
$D$
and length
$L$
injected with an impact velocity
$v_\perp =v_{\textrm {inj}}\sin \theta _{\textrm {s}}$
(component of the injection velocity
$v_{\textrm {inj}}$
normal to the shatter plane at an angle
$\theta _{\textrm {s}}$
), this fragmentation model yields the characteristic size distribution

Here,
$\alpha =X_{\textrm {R}}/D$
,
$\beta =X_{\textrm {R}}/LC$
and
$X_{\textrm {R}}=v_\perp ^2/v_{\textrm {thres}}^2$
, the latter of which is the ratio of the impact energy and the threshold energy required to initiate fracture in the pellet and
$K_0$
is the modified Bessel function of the second kind accounting for the energy propagation through the pellet. The material parameters
$C$
and
$v_{\textrm {thres}}$
depend on the neon–deuterium content ratio
$f_{\textrm {Ne}}$
of the pellet, both of which are obtained from Gebhart et al. Reference Gebhart, Baylor and Meitner(2020b). Fragments are sampled from (2.7) in sequence until the sum of all fragment volumes exceeds the pellet volume, after which the final fragment is shrunk such that the total volume equals that of the pellet. In doing so,
$N_{\textrm {frag}}$
will depend on the specific random seed used during the sampling process.
The next step concerns assigning velocities to each fragment. For simplicity, the magnitude and direction are sampled independently from one another and uncorrelated to the corresponding fragment’s size. The velocity magnitude
$v_{\textrm {frag}}$
is sampled from a normal distribution, where the mean fragment speed
$\overline {v}_{\textrm {frag}}$
is assumed to be given by the average of the injection speed
$v_{\textrm {inj}}$
and its component parallel to the shattering surface at an angle of
$\theta _{\textrm {s}}=12.5^\circ$
, i.e.
$\overline {v}_{\textrm {frag}}=v_{\textrm {inj}}(1+\cos \theta _{\textrm {s}})/2$
, and the spread is set to be
$\Delta v_{\textrm {frag}}/\overline {v}_{\textrm {frag}}=20\,\%$
(see Peherstorfer Reference Peherstorfer2022). Finally, each fragment is assigned a direction in which it traverses through the plasma at constant speed that is sampled uniformly within a circular cone with an opening angle of
$20^\circ$
(Peherstorfer Reference Peherstorfer2022).
2.4. Simulation set-up
Each simulation presented in this paper uses initial conditions from the representative AUG discharge no. 40655, which is an H-mode deuterium plasma with a major radius of the magnetic axis
${R_0}=1.74\,\textrm {m}$
, minor radius
$a=0.39\,\textrm {m}$
(midplane low-field side distance from
$R_0$
to the last closed flux surface) and on-axis toroidal magnetic field
$B_0=1.81\,\textrm {T}$
. Reconstructed magnetic equilibrium data and relevant initial plasma profiles are taken at
$t=2.3\,\textrm {s}$
, which is approximately the time of the SPI trigger, all of which are obtained from the integrated data analysis (IDA) framework (Fischer et al. Reference Fischer, Fuchs, Kurzan, Suttrop and Wolfrum2010, Reference Fischer2016, Reference Fischer, Giannone, Illerhaus, McCarthy and McDermott2020). The IDA framework applies Bayesian probability theory to the combined analysis of multiple diagnostics and physical modelling to obtain the most cohesive profile reconstructions and uncertainty estimates.
The reconstruction of the initial electron density profile
$n_e$
is based on both deuterium cyanide laser interferometry and Thomson scattering spectroscopy, and the electron temperature
$T_e$
is based on the latter. The initial thermal ion temperature
$T_i={T_{\textrm {D}}}$
is based on charge exchange recombination spectroscopy. Furthermore, neutral beam injection (NBI) provides a non-negligible population of super-thermal fast ions, for which IDA makes estimates of the fast ion density
$n_{\textrm {fi}}$
and pressure
$p_{\textrm {fi}}$
profiles based on kinetic modelling using the RABBIT code (Weiland et al. Reference Weiland, Bilato, Dux, Geiger, Lebschy, Felici, Fischer, Rittich and van Zeeland2018). With the fast ions being made up of ionised deuterium, and assuming the pre-SPI plasma is ideal, the quasi-neutrality criterion implies that the initial thermal ion density can be expressed as
$n_{\textrm {D}}=n_e-n_{\textrm {fi}}$
. The total plasma pressure
$n_eT_e+n_{\textrm {D}}{T_{\textrm {D}}}+p_{\textrm {fi}}$
is then used to solve for the magnetic equilibrium in the Grad–Shafranov equation, which is done iteratively since these profiles themselves depend on the equilibrium (see Fischer et al. Reference Fischer2016).
In this model we do not resolve the momentum space distribution of any of the particle species involved. Instead, the thermal and fast ions are incorporated into a single Maxwellian distribution, with density
$n_{\textrm {D}}=n_e$
and a rescaled temperature as to maintain the total plasma pressure. By doing this, it is effectively assumed that these fast ions instantaneously thermalise at the start of the simulation. It is also assumed that NBI is turned off at this point in time, meaning no further fast ions are introduced.
The corresponding flux surface geometry, with the SPI plume indicated, and the relevant plasma profiles are shown in figure 1. Note that the total current density from experiment includes bootstrap current and the NBI-driven current, which is here incorporated into the Ohmic current, as shown in figure 1(d). This is motivated by the fact that the bootstrap current quickly becomes Ohmic as a result of the relaxation of the pressure profile during the enhanced transport event. The NBI-driven current is turned off by the control system as soon as insufficient absorption is detected, and since we do not model the complex control system logic, this contribution to the total current is omitted.

Figure 1. Input data based on the reference AUG discharge no. 40655, used in all simulations presented. The illustration in (a) shows the flux surface geometry and the SPI plume, within which all fragments travel in straight lines originating from the shatter point. The particular shatter head considered in this study has a square cross-section with
$22\,\textrm {mm}$
side lengths, and a shatter angle
$\theta _{\textrm {s}}=12.5^\circ$
. Flux surfaces are indicated in red, and the
$q=2$
surface is marked in blue. The initial plasma profiles used in this work consist of (b) density and (c) temperature profiles of the bulk electrons (blue), thermal ions (red) and fast ions (green) and (d) the initial current density.
The radial coordinate is discretised on a grid with
$n_r=11$
points, and for the time coordinate we use an adaptive time step based on the assumption that the shortest time scales are those of ionisation
$\tau _{\textrm {ionis}}\sim |{\partial \log n_e/\partial t}|^{-1}$
. It sets the next time step according to
$\Delta t=\min (\Delta t_{\textrm {max}}, \alpha \tau _{\textrm {ionis}})$
, where we set
$\Delta t_{\textrm {max}}=10^{-6}\, \textrm {s}$
as the maximum allowed time step and
$\alpha =10^{-3}$
. With this set-up, depending on the pellet neon content, the typical wall-clock time for the simulations considered here is between 5 and 15 minutes on a computer running on an AMD EPYC 9354 processor.
2.5. Experimental cases selected for comparison
The SPI system installed for the 2022 AUG SPI experimental campaign is highly flexible, with three different shatter heads, of which two have a square cross-section at a
$\theta _{\textrm {s}}=12.5^\circ$
and
$25^\circ$
shatter angle, respectively, and the third has a circular cross-section at
$\theta _{\textrm {s}}=25^\circ$
. With two different shatter angles it is possible to isolate the normal impact velocity
$v_\perp =v_{\textrm {inj}}\sin \theta _{\textrm {s}}$
, which sets the fragment size distribution as given in (2.7), and the velocity distribution of the fragments
$v_{\textrm {frag}}$
. A detailed description of the SPI system is provided by Dibon et al. (Reference Dibon2023) and Heinrich et al. (Reference Heinrich, Papp, de Marné, Dibon, Jachmich, Lehnen, Peherstorfer and Vinyar2024).
In the 2022 AUG SPI campaign, approximately 180 SPI shots were carried out with this same H-mode scenario, of which no. 40655 was selected as the representative discharge. To represent a scan in injected neon quantity, we have chosen a small subset of discharges from the same scenario, each of which have relevant plasma parameters (such as plasma stored energy or pedestal top temperature) within
${\lesssim}5\,\%$
of the reference shot. In this paper, all the considered cases are pellets injected via the square shatter head with
$\theta _{\textrm {s}}=12.5^\circ$
(GT3). With input conditions all based on the reference AUG discharge no. 40655, as described in § 2.4, only the injection parameters are modified when modelling the various discharges, mainly the fraction of neon
$f_{\textrm {Ne}}$
in the injected pellet. An overview of the different set-ups is shown in Table 1. The values for the pellet length
$L$
, diameter
$D$
and injection speed
$v_{\textrm {inj}}$
are based on image processing analysis by Peherstorfer (Reference Peherstorfer2022). Note, that the 40 % neon case has a higher injection velocity than the rest. There were no experiments executed in 2022 with 40 % neon and a matching velocity, but we nevertheless wished to include a neon fraction case between 10 % and 100 % in the experimental comparison. For each discharge we run 40 otherwise identical simulations with varying random seeds, resulting in a variation in the number of fragments sampled between the simulations. The last column in Table 1 includes the average and standard deviation in the number of fragments.
Table 1. The SPI parameters used in the simulations of the considered AUG discharges, being the pellet neon fraction
$f_{\textrm {Ne}}$
, number of injected neon atoms
$N_{\textrm {Ne}}$
, pellet length
$L$
and diameter
$D$
and pellet injection speed
$v_{\textrm {inj}}$
. All injections use the
$\theta _{\textrm {s}}=12.5^\circ$
shatter head with a square cross-section – corresponding to guide tube 3 (GT3) in the 2022 experimental set-up – as indicated in figure 1(a). The final column shows the mean and standard deviation in the number of fragments generated using these parameters, based on 40 samplings of the fragment size distribution (2.7).

3. Simulation results and experimental comparison
3.1. Pellet neon fraction scan
To gauge the effect the pellet neon fraction
$f_{\textrm {Ne}}$
(or equivalently, the number of injected neon atoms,
$N_{\textrm {Ne}}$
) has on this dynamic system, we scan in the range between
$100\,\%$
down to
$0.0001\,\%$
in search of trends. Figure 2 shows the plasma current, electron temperature at the
$q=2$
surface and the total radiated power, as functions of time for different values of
$f_{\textrm {Ne}}$
. The initial temperature drop in all cases is due to dilution cooling as the density is rapidly increased. With more neon introduced to the system, more energy is lost via radiation, leading to a more rapid drop in temperature, which in turn causes an earlier onset time for the enhanced transport event triggered. Since the conductivity scales as
$\sigma \propto T_e^{3/2}$
, a stronger temperature drop with an increasing amount of injected neon accelerates the current drop during the current quench.

Figure 2. Scan in the amount of injected neon, showcasing the evolution of the (a) plasma current, (b) electron temperature at the resonant
$q=2$
surface and the (c) total radiated power. The dashed line in (b) indicates the temperature threshold of
$10\,\textrm {eV}$
at which we trigger the enhanced transport event. The scan consists of 300 individual simulations, each with fragments sampled from (2.7) which depends on
$f_{\textrm {Ne}}$
.
For some cases with less than
$5\,\%$
trace amounts of neon, or
$N_{\textrm {Ne}}\approx 1.5\times 10^{21}$
, rather than moving in from the edge towards the
$q=2$
surface, the thermal collapse below
$10\,\textrm {eV}$
is initiated in the core region. When the fragments pass
$R={R_0}$
and reach the high-field side, they again deposit material on flux surfaces they already passed on the low-field side. For cases with neon fractions
$f_{\textrm {Ne}}\lesssim 0.001\,\%$
, or
$N_{\textrm {Ne}}\lesssim 3\times 10^{17}$
, a thermal collapse does not occur, and the plasma therefore does not disrupt.
The reason for the deep fragment penetration in the simulations is twofold. Less neon in the simulations will allow the fragments to penetrate deeper because it takes more time before the thermal collapse occurs due to less radiative cooling. At the same time, with smaller neon fractions the Parks model predicts fewer fragments that take longer to ablate, further improving their penetration. In this manuscript we do not isolate the effects of pellet composition from the number of fragments. A similar scan in the pellet composition is done by Patel et al. (Reference Patel2025), in which the number of fragments is kept fixed at 200 across all neon fractions. This same phenomenon of an inside-out collapse at lower neon fractions is also observed, suggesting it is in fact the lower radiative cooling that allows for a deeper fragment penetration, and not only the low number of fragments generated. Core penetration of SPI is also observed experimentally in the ultra-high-speed video recordings of AUG SPI experiments (Patel Reference Patel2023).
In the context of runaway electron generation, the maximum simulated runaway current is approximately
$5\,\textrm {A}$
, occurring in the cases with majority neon admixture. This is mainly attributed to the larger runaway seed generated by the hot-tail mechanism due to the more pronounced thermal quench. However, the induced electric field following this drop in electron temperature is insufficient to produce a significant amount of runaway electrons from this seed population via the avalanche mechanism.
3.2. Plasma current evolution
Simulated plasma currents of six selected injection scenarios are compared with experimental measurements in figure 3. Each case was run 40 times using different realisations of the fragment sampling procedure. For the cases with a finite amount of injected neon, the simulations generally agree well with experiment. Note that since the magnetic field is kept static in the simulations, no runaway scrap-off is accounted for (Vallhagen et al. Reference Vallhagen, Hanebring, Fülöp, Hoppe and Pusztai2024b ), and effects of vertical displacement events (VDEs) are not captured. These VDEs may explain the kink seen in the experimental data at approximately half-way down the current quench. In addition, we do not model any control system action, as the AUG discharge control system is too complex to implement in such a simulation, and a feed-forward implementation of the actions executed would lead to a departure from self-consistent modelling of the plasma evolution.

Figure 3. Comparison of the plasma current evolution for varying amounts of injected neon in experimentally measured data (black) and simulations (red), for which each case is run 40 times using different random seeds while sampling fragments. Except for the pure deuterium injection, a good agreement with experimental plasma currents is observed during the initial current quench evolution (100 %–80 % of
$I_{\textrm {p}}$
) prior to the onset of the VDE phase.
The last column in Table 1 displays the average number and standard deviation of fragments sampled for the six injection scenarios. The fragment size distribution (2.7) has a dependence on the neon fraction and the impact velocity
$v_\perp$
(or
$v_{\textrm {inj}}$
, as we fix the shatter angle
$\theta _{\textrm {s}}=12.5^\circ$
) that determines the size of the fragments, or equivalently the total number of fragments sampled.
In discharge no. 40738 no neon is introduced into the plasma, and in the simulations the radiated losses are not sufficient to initialise the enhanced transport event, meaning that the simulated plasma does not disrupt. This was expected in the modelling due to the relatively low radiation of deuterium. In the experiment we observed a relatively long quench of the plasma current. Full-sized pure deuterium pellets injected into the AUG SPI scenario plasma typically lead to a disruption, either via an inadvertent control system action (Sieglin et al. Reference Sieglin, Maraschek, Kudlacek, Gude, Treutterer, Kölbl and Lenz2020; Heinrich et al. Reference Heinrich2025), or, as is suspected in this case, due to background impurity radiation (mainly tungsten from the first wall). The presence of background impurities provides additional loss channels through line radiation, enabling the decrease in temperature below
$10\,\textrm {eV}$
within the
$q=2$
flux surface, and thus triggering the enhanced transport event. However, it was also observed that 100 % D
$_2$
injections containing fewer than
$N_{\mathrm{D_2}}=10^{22}$
deuterium molecules do not disrupt the H-mode scenario used for the SPI experiments (Heinrich et al. Reference Heinrich2025).
3.3. Radiated energy fraction
The radiated energy fraction
$f_{\textrm {rad}}$
shows how much of the total available plasma energy drop is lost during a disruption via radiation, as opposed to conduction. It is calculated using the formula by Lehnen et al. (Reference Lehnen2011) and Sheikh et al. (Reference Sheikh2020)

where
$W_{\textrm {rad}}$
is the radiated energy,
$W_{\textrm {th}}$
and
$W_{\textrm {mag}}$
are respectively the initial thermal and poloidal magnetic energy,
$W_{\textrm {ext}}$
is external heating and
$W_{\textrm {c}}$
is the energy coupled into the surrounding conducting structure. For the simulations, the radiated energy is calculated as
$W_{\textrm {rad}}=\int _{t_{\textrm {start}}}^{t_{\textrm {end}}}\textrm {d}{t}P_{\textrm {rad}}$
, where the time
$t_{\textrm {start}}$
is when the first fragment reaches the plasma, and
$t_{\textrm {end}}$
is taken to be at the end of the radiation peak, when the derivative of
$P_{\textrm {rad}}$
approaches zero, as was done for the experimental estimates (Heinrich et al. Reference Heinrich2025). The initial thermal energy
$W_{\textrm {th}}$
is the combined energy of all thermal particles, i.e.
$W_{\textrm {th}}=3/2(n_eT_e+\sum _i{n_i}T_i)$
, evaluated at
$t=t_{\textrm {start}}$
. Typically in AUG experiments, the heating systems automatically turn off due to insufficient absorption shortly after the fragments enter the plasma, and will therefore contribute to the total energy of the system for a brief time. In the experimental evaluation of
$f_{\textrm {rad}}$
we correct for this by subtracting
$W_{\textrm {ext}}$
, which accounts for the external heating of the plasma between the arrival of fragments and shut-off of the external heating systems (Heinrich et al. Reference Heinrich2025). In the current simulation set-up, we can simply assume that all heating system are effectively turned off at
$t=t_{\textrm {start}}$
, and thus we set
$W_{\textrm {ext}}=0$
. Furthermore, the energy coupled to the surrounding structure
$W_{\textrm {c}}=0$
, since we assume a perfectly conducting wall as described in § 2.1.
Figure 4 shows a comparison between the radiated energy fraction from simulations and experimental estimates, as described by Heinrich et al. (Reference Heinrich2025). The calculation involves interpolating foil bolometry measurements of the radiated power in different sectors of the tokamak, and integrating it over an interval of time similarly determined as in the simulation calculations. The experimental values have an estimated relative uncertainty of approximately
$20\,\%$
.

Figure 4. Radiated energy fraction
$f_{\textrm {rad}}$
from simulations (red) and experimental estimates (black) as functions of the amount of injected neon. The scatter in the simulation data comes from generating fragment samples using 40 different random seeds in otherwise identical simulations. For the experimental estimates, a relative uncertainty of
$20\,\%$
(grey) is assumed (Heinrich et al. Reference Heinrich2025).
Another key difference between simulation and experiment for 100 % D
$_2$
injection, apart from the plasma current, is
$f_{\textrm {rad}}$
. Otherwise, the simulated radiated energy fraction aligns well with experiment results, both in magnitude and its scaling with the neon fraction. While within the experimental uncertainty margins,
$f_{\textrm {rad}}$
for the highest neon fractions is close to 100 % as opposed to the 80 % observed in the experiments. There is yet no definite conclusion for the reason of the saturation of the experimentally observed
$f_{\textrm {rad}}$
around 80 %. A potential cause for this may be the incorrect estimation of conducted and coupled energy losses. However, both of these are difficult to self-consistently model in a one-dimensional simulation.
Compared with three-dimensional MHD simulations (Tang et al. Reference Tang2025) we see that the one-dimensional Dream simulations somewhat overestimate
$f_{\textrm {rad}}$
. While Dream provides a better match than Jorek for the experiments at low neon doping above approximately 1 %, neon content Jorek simulations show a radiated energy fraction closer to the experimental values than Dream.
4. Discussion and conclusion
We applied a one-dimensional fluid model of SPI-induced disruptions in the ASDEX Upgrade tokamak, dynamically evolving the plasma current, thermal bulk of electrons and ion charge state densities. This numerically inexpensive model was demonstrated to yield good results when comparing with experimental data, capturing the general trend in the plasma current evolution and radiated energy fraction when varying the amount of injected neon. We also observe no significant runaway generation, as in the experiments.
A notable exception from the close agreement is the case of 100 % deuterium injection, where the experiments do disrupt eventually (if
$N_{\mathrm{D_2}}\gtrsim 10^{22}$
), typically due to the combined effect of background impurity radiation, reaching the density limit, vertical instability or control system action. These aspects are not yet modelled here. In particular, the discrepancy between the simulated radiated energy fraction and the experimental estimates is primarily attributed to the absence of background impurities in the simulations. These impurities would enhance radiative energy losses, reducing the temperature sufficiently to trigger the observed enhanced transport event. The inclusion of tungsten as a background impurity in the plasma as well as the material deposition shift due to plasmoid drifts is proposed for future investigations.
A probabilistic method for generating individual fragment sizes and velocities is presented, the former of which is based on fragmentation model first proposed by Parks (Reference Parks2016) and modified by Gebhart et al. (Reference Gebhart, Baylor and Meitner2020a ,Reference Gebhart, Baylor and Meitner b ). With the present set-up based on Parks, statistical variations in this sampling procedure provide an uncertainty quantification and are shown to only have a notable impact on the disruption dynamics for trace amounts of injected neon. The same is also observed experimentally (Heinrich et al. Reference Heinrich2025). This means that, while thermal load mitigation would be difficult to optimise by selecting appropriate size distributions, it is also not a factor to worry about when determining the optimal SPI parameters for e.g. material assimilation for runaway electron suppression.
A potential cause for discrepancies between simulation and experiment is that sampling from the Parks fragment size distribution has been shown by Peherstorfer (Reference Peherstorfer2022) to differ from the experimentally measured fragment size distributions at AUG. The Peherstorfer analysis was, however, limited to
$4\,\textrm {mm}$
diameter pellets with typically low neon content, most of which do not cause disruptions at AUG, and thus are not suitable for experimental validation of disruption simulations. Furthermore, we wanted to first compare results obtained with the Parks model with the experiments, as this is widely used in the literature and thus a validation exercise was considered useful for the community. Recent peridynamic modelling by Lee et al. (Reference Lee, Madenci, Na, de Marné, Dibon, Heinrich, Jachmich, Papp and Peherstorfer2024) has shown a promising match between modelled and measured fragment size distributions. Furthermore, reconstructed fragment size and velocity distributions from fast camera videos using a deep learning-based algorithm (Illerhaus et al. Reference Illerhaus, Treutterer, Heinrich, Miah, Papp, Peherstorfer, Sieglin, Toussaint, Zohm and Jenko2024) are becoming available. Investigating the effect of these different fragment distributions on the disruption dynamics is planned for a future study. Comparisons with three-dimensional MHD simulations can help to inform model parameters and joint studies are planned in the future with fast one-dimensional scans complemented by more numerically expensive nonlinear three-dimensional MHD studies of the most relevant scenarios (Tang et al. Reference Tang2025).
Acknowledgements
The authors are grateful to S. Jachmich, J. Artola, W. Tang and A. H. Patel for the fruitful discussions.
Editor Tünde Fülöp thanks the referees for their advice in evaluating this article.
Funding
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). 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 the ITER Organisation. Neither the European Union nor the European Commission can be held responsible for them. The ASDEX-Upgrade SPI project has been implemented as part of the ITER DMS Task Force programme. The SPI system and related diagnostics have received funding from the ITER Organisation under contracts IO/20/CT/43-2084, IO/20/CT/43-2115, IO/20/CT/43-2116.
Declaration of interests
The authors report no conflict of interest.