Modelling of runaway electron dynamics during argon-induced disruptions in ASDEX Upgrade and JET

Disruptions in tokamak plasmas may lead to the generation of runaway electrons that have the potential to damage plasma-facing components. Improved understanding of the runaway generation process requires interpretative modelling of experiments. In this work we simulate eight discharges in the ASDEX Upgrade and JET tokamaks, where argon gas was injected to trigger the disruption. We use a fluid modelling framework with the capability to model the generation of runaway electrons through the hot-tail, Dreicer and avalanche mechanisms, as well as runaway electron losses. Using experimentally based initial values of plasma current and electron temperature and density, we can reproduce the plasma current evolution using realistic assumptions about temperature evolution and assimilation of the injected argon in the plasma. The assumptions and results are similar for the modelled discharges in ASDEX Upgrade and JET, indicating that the implemented models are applicable to machines of varying size, which is important for the modelling of future, larger machines. For the modelled discharges in ASDEX Upgrade, where the initial temperature was comparatively high, we had to assume that a large fraction of the hot-tail runaway electrons were lost in order to reproduce the measured current evolution.


Introduction
Runaway electrons (RE) may cause severe damage to plasma-facing components in tokamaks [1], where they may occur as a consequence of disruptions. To avoid this, different schemes are being developed to mitigate, limit or entirely avoid the formation of REs. Massive material injections (MMI) are proposed to this end, which may be realized through a gas injection (massive gas injection -MGI) or the injection of solid pellets, which can be shattered when entering the vacuum chamber (shattered pellet injection -SPI) [2,3]. In medium-sized tokamaks, the potential of these measures have been demonstrated [4][5][6][7][8][9][10][11][12][13][14], but in the future devices envisaged for fusion energy generation, the plasma conditions will be significantly different (higher temperatures and densities, larger plasma currents), and whether the proposed measures for RE mitigation or avoidance will be effective also in these devices can currently only be assessed through plasma physics modelling.
Several theoretical models for the physics of runaway electrons in tokamak disruptions have been developed and implemented in computational tools [15]. To assess to which extent these models can be applied to make useful predictions for the plasma behaviour during disruptions, they must be validated against existing experimental data.
In current tokamaks, disruptions are often deliberately triggered by MGI, leading to a thermal quench (TQ) which under certain conditions is followed by a formation of REs and a rapid decay of the ohmic plasma current (a current quench -CQ). The experimental data collected during such discharges constitutes a valuable dataset for validation of models for the formation of REs in the presence of impurities (often in the form of massive amounts of noble gasses). In this work, we use such experimental data from the tokamaks ASDEX Upgrade (AUG) and Joint European Torus (JET) to investigate the applicability of the models, possible modifications needed, and qualitative differences between the mechanisms behind RE formation in these two differently sized tokamaks.
The main computational tool used in the present paper is called go, a fluid code that describes the radial dynamics of the current density and the electric field, in the presence of impurities. The models implemented in this tool are briefly described in section 2, and more details are given in Refs. [16][17][18][19]. In the version of go that is used in the paper, hot-tail, Dreicer and avalanche RE generation models are implemented, as will be further described in section 2.
REs are defined as electrons having a momentum larger than the critical momentum p c (or, equivalently, a velocity larger than the critical velocity v c ), above which the accelerating force exerted by the electric field (induced in the disruption) is larger than the friction force due to collisions, i.e. electrons above this threshold are accelerated until their momentum is limited by radiative energy losses. When the collision time at the critical momentum is longer than the duration of the TQ, hot-tail is the dominant primary generation mechanism [20]. In future fusion devices, this is expected to be the case [21], so significant hottail generation is expected in e.g. ITER [22] or SPARC [23]. For this reason, it is important to properly understand and being able to model the interplay of the various runaway generation mechanisms in experimental scenarios.
Kinetic simulations of AUG discharges, using the tool code [24], were recently presented by Insulander Björk et al. [25]. code solves the linearized Fokker-Planck equation including radiation reactions and avalanche source, as well as the electric field evolution [26]. Kinetic simulations of JET discharges with code were attempted, but became prohibitively slow due to the high electric fields induced. Also, code lacks radial diffusion of electric field and current, mechanisms which have been shown to be important for modelling of disruptions in large machines such as ITER [19], so in this sense, the two codes are complementary.
In this paper, we use the go code to investigate argon-induced disruptions in the AUG and JET tokamaks. In agreement with recent results based on coupled fluid and kinetic simulations with go and code [27], we find that taking into account all the hot-tail electrons in AUG would overestimate the final runaway current. To be able to match the experimentally obtained current evolution, we vary the loss-fraction of the hot-tail RE. Such losses are expected to occur due to the breakup of magnetic flux surfaces, accompanying the TQ.
Furthermore, we find that the self-consistent calculation of the temperatures of ions and electrons using time-dependent energy transport equations often results in a temperature evolution which does not agree with measurement data. The main reason why the energy transport equations fall short of accurately predicting the temperature evolution is likely the lack of a detailed self-consistent model for the threedimensional evolution of the magnetic field during the TQ. However, by assuming an exponentially decaying temperature evolution throughout the simulation, the experimentally observed current evolution can be reproduced both in AUG and JET.

Fluid modelling with go
go [16][17][18][19] simulates the radial dynamics of the of the current, temperature, density and electric field. Models for the previously mentioned RE generation mechanisms are implemented in the code, as well as RE generation through Compton scattering and tritium decay, albeit the two latter mechanisms are not relevant for the simulations presented here.

Runaway electron generation
The hot-tail RE generation is a result of a rapid cooling of the plasma, during which there is not enough time for the fastest electrons (the "hot tail" of the electron velocity distribution) to thermalize, thereby ending up above the critical momentum. The analytical model used in go for hot-tail generation is derived in Ref. [28].
t 0 n(t)/n 0 dt, n is the electron density, v T is the thermal electron speed, v c is the critical velocity, subscript 0 denotes an initial value, H is the Heaviside function, ν = ne 4 ln Λ/(4πǫ 2 0 m 2 e v 3 T ), lnΛ is the Coulomb logarithm, ǫ 0 is the vacuum permittivity and m e the electron rest mass. This expression takes into account the directivity of the electric field, in contrast to the simplified expression used in Refs. [22,29], which is based on counting the number of electrons with v > v c .
The Dreicer RE generation is instead a consequence of momentum space diffusion of electrons to momenta above the critical threshold, due to small angle collisions [30]. For modelling of Dreicer RE generation, a neural network was used, trained on the output of kinetic simulations of plasmas with impurities, which were performed with code [31].
RE generation through these two mechanisms is referred to as primary generation. Already existing REs thus generated may transfer part of their momentum to bulk electrons, knocking them above the critical momentum and thereby creating new REs at an exponentially growing rate. This a secondary generation mechanism, referred to as avalanche [32,33]. The avalanche generation is modelled by a semi-analytical formula [34], that has been carefully benchmarked against kinetic code simulations of disruptions in impurity-containing plasmas.

Temperatures
The temperature evolution in this work was simulated using different models, described below.
2.2.1. Exponential temperature model During the initial phase of the simulations when the temperature decays rapidly (TQ), the magnetic flux surfaces are broken up due to magneto-hydrodynamic (MHD) instabilities [35]. The rapid radial transport due to these effects is the main driving force behind the initial temperature change. This complex transport process cannot be self-consistently modeled within our fluid framework, and emulating it in terms of spatiotemporally varying transport coefficients would introduce a large set of unconstrained free parameters. Instead, in our simulations the on-axis temperature T oa is initially prescribed to follow an exponential decay [22,28]: where the initial temperature T initial is deduced from experimental data, and the thermal quench time t TQ and the temperature at the end of this phase and most of the CQ, T CQ , may be treated as free parameters which can be adjusted to yield an I p evolution similar to the measured data. Another reason for using this prescribed temperature evolution during the TQ is to be consistent with the assumptions used to derive the analytical model for the hot-tail runaway generation [28]. During this initial temperature drop, the radial temperature profile shape is assumed to remain constant in the present simulations. The initial radial temperature profile T initial is deduced from experimental data as described in sections 3.1.1 and 3.2.1, for AUG and JET respectively, and this distribution is then scaled with the value given by equation (2) to calculate the free electron temperature in each radial point. Whereas this is undoubtedly a coarse approximation, experimental measurements cannot yield useful data on the actual temperature profile during this brief and turbulent phase of the disruption to support a better estimate. Assuming a flat temperature profile with a value corresponding to a homogeneous spreading of the initial thermal energy will lead to slightly higher final runaway currents, but do not lead to qualitative changes in the results. Note, that since T CQ is assumed to be constant, the final temperature profile is flat in all simulations.
2.2.2. Energy transport temperature model go also has the capability to calculate the temperatures of free electrons, bulk ions and impurity ions using time-dependent equations for the energy transport of the respective species in a one-dimensional cylindrical model [17]. Whereas these equations do not account for particle transport, and are not sufficiently well constrained during the TQ, an attempt was made to use them for calculation of the post-TQ temperatures.
The energy transport equations are implemented as where n is the total density of all species (electrons and ions), σ is the conductivity, Z eff is the effective ion charge, n k i is the density of the i th charge state (i = 0, 1, ..., Z − 1) of the ion species k (e.g. deuterium, argon), P Br and P ion are the Bremsstrahlung and the ionization energy losses, respectively. The line radiation rates L i k (T, n e ) are extracted from ADAS [36]. In these calculations, which were only attempted for the post-TQ plasma, a constant heat diffusion coefficient χ = 1 m 2 /s was used, and it was confirmed that the results were insensitive to this choice.
In these simulations, the temperature is assumed to be equal for all species, since this is computationally less expensive. This simplification affects the heat capacity with a factor of at most two [19]. Since the initial part of the temperature evolution, where most of the thermal energy is lost, is not simulated with this model, this simplification is not expected to have an important impact. The injected argon is instead assumed to be heated through heat exchange with the particles present before the injection.

2.2.
3. Hybrid temperature model An attempt was made to use a hybrid temperature model. First, the temperature was calculated using equation (2) with T CQ = 1 eV, until T oa reached a pre-set value T switch .
It has been estimated [37] that transport due to stochastic flux surfaces could drive a temperature drop down to ∼100 eV more effectively than radiation, which provides a physically motivated value for T switch . Then the simulation was restarted using the time-dependent energy transport equation (3) for temperature determination, i.e. the value T CQ = 1 eV is not reached. After the switch to using the energy transport equations the temperature is calculated locally in each radial point by go.
In many cases, this strategy resulted in a rapid increase in the calculated on-axis temperature up to a few 100 eV after the switch. Re-heating has been sporadically observed in natural disruptions [38], but has not been possible to investigate experimentally in MGI-induced disruptions. On ASDEX Upgrade, the electron cyclotron emission temperature diagnostic (see section 3.1.1) is in density cut-off after the MGI, therefore any potential re-heating can not be directly measured. However, it is likely that the losses in the energy transport equation are underestimated, for example due to radiation from wall impurities or remaining transport losses that are not included, or underestimated, in the present model. The radiation losses with argon impurities present follow a nonmonotonic behaviour as a function of temperature, with maxima at ∼10 eV and ∼100 eV. Even a small underestimate of the radiative losses can make the ohmic heating overcome the cooling in the ∼10 eV range, and make the temperature evolve towards the ∼100 eV range where the ohmic heating is less efficient. However, if the calculated losses are large enough, the temperature will instead move towards the ∼10 eV range, where it remains until the ohmic heating has decayed. When the simulations result in such a temperature evolution also the experimental I p evolution was reproduced fairly well.
There is a close connection between the current decay rate and the temperature through the conductivity. In cases where re-heating occurred in the simulations, the experimentally observed current decay rate was not reproduced sufficiently well. Thus, in sections 4.2 and 4.4 we used the exponential temperature model only, i.e. equation (2) was used for the entire simulation, and the temperature T CQ was adjusted, together with other free parameters, to match the measured I p evolution. The values of T CQ which could match the I p evolution were approximately 20 eV, i.e. T CQ is the free electron temperature during the most of the CQ, but not necessarily the final temperature during the RE plateau phase.

Densities
go calculates the ionization states of the argon and the resulting free electron density from the time-dependent rate equations using as input the free electron density prior to the injection, the injected amount of argon and the initial density profile. I i k denotes the electron impact ionization rate and R i k the radiative recombination rate for the i th charge state of species k, respectively. The ionization and recombination rates are extracted from ADAS [36]. The injected argon is assumed to distribute according to the same density profile as the initial deuterium density.
In the current simulations, the quantity of assimilated argon is given in terms of the ratio r Ar/D of the argon density n Ar to the initial deuterium density, assumed to be equal to the initial free electron density n e0 . Lacking reliable experimental data on the rate of assimilation of argon into the plasma, the argon is assumed to have assimilated and distributed within the plasma before the beginning of the simulated time span. The density of free electrons and of each ionization state is calculated by go. T oa , calculated by equation 2, is the on-axis free electron temperature including cooling by the injected argon. The effect of modelling the argon density as exponentially increasing after the beginning of the simulation was investigated, but this introduced another free parameter (the time rate of argon assimilation) which affected the simulation results in a similar manner as assuming a constant r Ar/D , so the simpler approximation was preferred.
The evolution of the free electron density during the disruption is difficult to diagnose directly. For AUG, the total assimilated amount of argon was deduced from the current dissipation rate during the RE plateau phase, as explained in Ref. [25], and the same assimilation rates were used in the present simulations. For JET, experimental data was used to determine the total amount of assimilated argon, as explained in section 4.4. However, we did not attempt to deduce the time evolution of the argon assimilation from these data, partly due to their uncertain nature, and partly in order to use the same modelling strategy for both tokamaks.

Plasma current
go models the evolution of the current density as a balance between the diffusion of the electric field and the generation of the REs: where j = σ Sp E + ecn RE is the sum of Ohmic and runaway current densities, with σ Sp being Spitzer conductivity, and µ 0 , e, and c denote the magnetic permeability, the elementary charge and the speed of light, respectively. Toroidicity effects are neglected here.
In the disruptions modelled in this paper, the measured plasma current I p typically displays a small peak coinciding with the beginning of the TQ, which indicates an MHD event during which the current density is redistributed radially due to a breakup of magnetic flux surfaces [35]. Towards the end of the TQ, the magnetic flux surfaces re-form and the current can be confined until it is dissipated resistively.
The rapid relaxation of the current profile due to non-axisymmetric MHD is not modelled by go, and thus the calculated I p does not display this peak. After the peak, I p starts to decay rapidly (the CQ) due to the increase in plasma resistivity associated with the decreasing temperature (σ ∝ T 3/2 e ), which is modelled by go. Due to this temperature dependence, the I p evolution is strongly affected by T e evolution, and hence by the parameters describing it: T initial , T CQ and t TQ . As previously stated, the exact value of T CQ does not affect the results significantly if we use the hybrid temperature model, i.e. if we switch to the time-dependent energy transport equation before T CQ is reached. However, the value becomes important when we use the exponential decay temperature model, i.e. if equation (2) is used throughout the simulation, as discussed in sections 4.3 and 4.2. In these cases, T CQ rather plays the role of the equilibrium temperature that marks the end of the TQ, but not necessarily the final temperature reached after the CQ when no ohmic current remains causing resistive heating of the plasma.

Plasma elongation
go has the ability to model elongated plasmas [39], however, this capacity was not used in the presented simulations.
The modelled discharges featured slightly elongated plasmas, κ ≈ 1.16 in the JET discharges, according to pre-disruption equilibrium reconstructions ¶.
Inclusion of the elongation in the simulations did not result in any qualitative changes in the simulation results, and thus it was decided to reduce the parameter space by omitting this parameter. This simple geometry also facilitates comparison with previous kinetic simulations which were zero-dimensional in real space.

Experimental data
In order to assess the applicability of the go model for tokamaks with different parameter ranges, we model discharges in AUG and JET, and compare our results to experimental data. Table 1 shows some basic data typical for the discharges simulated for the respective tokamaks, and in the following sections, more specific data for each simulated discharge are presented.

ASDEX Upgrade
ASDEX Upgrade is a medium sized tokamak, and we model four discharges which were specifically designed for the study of runaway electron dynamics [5,12]. The runaway discharges studied in this paper are nearcircular (elongation κ ≈ 1.1), inner wall limited, have core ECRH (Electron Cyclotron Resonance Heating) and low pre-disruption density (n e0 ≈ 3·10 19 m −3 ). The discharges are terminated using argon MGI from an in-vessel valve. A 30 ms time period, starting ¶ EFIT/ELON at the argon valve trigger, was simulated to ensure that the entire current quench was covered in all the simulated cases. The plasma position remained radially and vertically stable in the experiments during the simulated time window. Four discharges with different plasma parameters and amounts of injected argon were selected for modelling and an overview of the basic parameters for these discharges is presented in table 2. These discharges were selected from the set of 11 discharges modelled earlier in Ref. [25].

Temperature
The temperature and its radial distribution in the plasma is measured using electron cyclotron emission (ECE). For the present simulations, we only use the initial temperature profile. The initial on-axis temperature T initial used in equation (2) is determined as the average temperature over the circular area within 1 cm of the magnetic axis and the time interval between 1 ms and 1.5 ms after the argon injection valve trigger. The time interval was selected to exclude both the beginning of the TQ and any initial temperature variations due to the shut-off of the heating system shortly before the disruption.

Densities
The free electron density is measured by CO 2 interferometry, which yields the line integrated free electron density along the line of sight of the interferometer. The density profile is fitted by the tool augped, which fits a modified hyperbolic tangent function [40] to the radial density data points obtained with interferometry.
The average measured free electron density during the first 1.5 ms after the argon valve trigger is used as the initial on-axis free electron density n e0 , since after these 1.5 ms, the argon injection causes the measured density to start to increase. The initial density profile is scaled with n e0 to obtain the initial density in each radial data point.
In each discharge, argon gas was injected into the vacuum vessel to trigger the disruption [5,41]. The injected number of argon atoms N Ar is calculated from the pressure in the MGI reservoir holding the argon gas before injection, the reservoir volume (0.1 l) and the gas temperature (300 K). This quantity is listed in table 2 for the respective discharges. The amount of argon which finally assimilates in the plasma as a function of time is difficult to assess. In previous work [25], the Ar assimilation fraction f Ar was assumed to be the same for all discharges, and was assessed by matching the current decay rate during the RE plateau phase. It was found that a reasonable estimate was that by the end of the TQ, 20% of the injected argon was assimilated in the plasma volume V p , as defined by the minor and major radii listed in table 1. We hence make the same assumption for the present simulations, and calculate the argon density n Ar as N Ar /V p · f Ar = N Ar /V p · 0.2. The ratio r Ar/D is calculated as n Ar /n e0 .

Joint European Torus -JET
JET is significantly larger than ASDEX Upgrade, as shown in table 1. Four discharges in which runaway electrons were formed were chosen for modelling, representing a set of varying plasma parameters. In particular, the size of the Ar injection varies by more than an order of magnitude. The basic parameters for these discharges are presented in table 3.

Temperature
The free electron temperature is measured by high resolution Thomson scattering (HRTS) which yields both a temperature profile and an sufficiently reliable value of the on-axis temperature before the Ar injection (T initial used in equation (2)). The initial temperature profile is smoothed using the rloess algorithm in matlab [42] to remove signal noise.

Densities
The free electron density n e is measured by several different diagnostics. Interferometry yields the line integrated free electron density along the line of sight of the interferometer. In these simulations, we used the fast interferometer signal + . The shape of + KG4C/LDE3 signal, compensated for fringe jumps and scaled to give a pre-disruption density in agreement with the DF/G1C- Table 2: Basic parameters for the four simulated discharges in AUG, and the notation for these parameters used in this paper. Further descriptions of the origins of these parameters are found in sections 3.  Table 3: Basic parameters for the four simulated discharges in JET and the notation for these parameters used in this paper. * In discharges #92454 and #92460, the Ar injection was succeeded by a large Kr injection. However, the effect on the I p evolution is insignificant until a few ms after the end of the modeled CQ. the initial free electron density profile was retrieved from the HRTS data at the last available data point for this diagnostic before the Ar injection, which was at most 130 ms before the injection. The profile was smoothed to remove signal noise using rloess [42] just as for the temperature profile, and then scaled so that the integral of the profile yielded the same initial line integrated density as the interferometry based density data. The final number of injected argon atoms is calculated from the volume of the reservoir holding the argon before the injection (DMV3, with volume 0.35 L), the gas vessel pressure before and after the injection * , and the assumption that the gas is held at room temperature (300 K). The fraction of this argon which actually assimilates in the plasma is difficult to assess experimentally, and the gas also leaves the injection reservoir during a period of some tens of ms. Hence, the ratio r Ar/D was regarded as a free parameter, constrained by the condition that the calculated maximum line integrated free electron density during the simulated time span should not deviate by more than 10% from the maximum line integrated free electron density given by the interferometry.

AUG simulations with the hybrid temperature model
In the first set of simulations, four AUG discharges were simulated using the same modelling parameters t TQ and r Ar/D as earlier [25] to verify that the modelled quantities show a qualitatively similar behaviour for go and the kinetic tool (code) used in Ref. [25]. For these simulations, the hybrid temperature model described in section 2.2.3 was used, i.e. equation (2) was used initially, and then we switched to the time-dependent energy transport equation after the on-axis T e had dropped to below a pre-defined temperature T switch . This method most closely resembles the strategy used in Ref. [25]. We set T switch = 6 eV, which is similar to the temperature equilibrium value found in several of the simulations in Ref. [25]. Furthermore, we set T CQ = 1 eV, and note that the value of T CQ has a minor impact on the exponential function, as long as it is lower than T switch , since the exponential function is abandoned before T CQ is reached. Relevant parameters from the simulations are shown in table 4 for all four discharges, along with the measured parameters, where applicable. The CQ time listed in table 4 is defined as the time from the beginning of the I p spike to the point where dI p /dt < 25 kA/ms, and the post-CQ total plasma current as the measured or calculated value of I p at this time point.
The calculated currents and RE generation rates are shown in figure 1 for AUG discharge #34183, as an example.
As explained in section 2.4, the phenomena resulting in the measured current spike are not modelled by go, and hence, the current spike is not seen in the plot of the simulated I p evolution. In addition to this expected discrepancy, we observe that in all the simulated discharges, the simulations • predict a partial CQ (i.e. neither full conversion nor complete CQ), • overestimate the total RE generation, leading to an underestimation of the CQ time and an overestimation of the post-CQ I p • overestimate the I p decay rate during the CQ and • predict a tiny Dreicer RE seed generation, relative to the hot-tail RE seed.
These qualitative indicators of the model's performance are identical to those reported in Ref. [25], i.e. we must conclude that these models need to be amended to be able to reproduce the observed I p evolution. The temperature evolution and the absence of RE losses are the main common features of these otherwise very different simulations, and in the following, we try to improve the agreement with experiment by changing to the exponential temperature model, and by modelling RE losses where appropriate.

AUG simulations with the exponential temperature model
As discussed in Ref. [25], the (almost) consequent overestimation of the post-CQ total plasma current is likely due the fact that losses of the RE seed population caused by the stochastization of the magnetic flux surfaces during the TQ are not modelled by code (and also not in go). In the code simulations reported in Ref. [25], the 2D momentum distribution of the entire electron population is modelled as a continuum, i.e. no explicit distinction is made between bulk electrons and RE. In contrast, go models the REs as a distinct electron population, and also models each RE generation mechanism separately. This makes it possible to model the loss of REs generated by the hottail mechanism by multiplying the number of hot-tail REs generated in each time step by a damping factor f d . For clarity, we discuss the hot-tail seed survival fraction f HT = 1 − f d in the following. The loss of hottail generated REs is expected to be significant since they are generated during the TQ when confinement is impaired due to the magnetic flux surfaces being broken up as noted in section 2.2.1. The losses of REs generated through the Dreicer and avalanche mechanisms are expected to be less important since (a) (b) Figure 1: (a) Current evolution and (b) RE generation rates for AUG discharge #34183, using the hybrid temperature model. In panel (a), the experimentally measured total plasma current is plotted for reference. The vertical lines mark the end of the CQ for the measured and calculated I p evolution respectively, see table 4. Note that in panel (b), the Dreicer generation rate is scaled by a factor 10 11 to be visible.
they are predominantly generated during the later phase of the disruption when the flux surfaces have re-formed. The four AUG discharges were simulated with hottail losses and the hybrid temperature model, but the overprediction of the post-CQ I p remained also when assuming loss of all hot-tail RE. When the hot-tail seed was assumed lost, the Dreicer generation increased instead, still resulting in a high total RE generation. In these cases, it was also found that the I p evolution depended strongly on the temperature at which the switch was made from using equation (2) to invoking the energy transport equation, while the model is only predictive of the temperature evolution during the CQ if it is insensitive to T switch in the physically motivated 10-100 eV range. Moreover, the calculated temperature in many cases increased to about 100 eV after the switch, which contradicts measurement data, as discussed in section 2.2.3. We therefore turned to the exponential temperature model, using equation (2) for the entire simulation, although other RE loss channels may also be part of the explanation for the inability to match observations. As mentioned in section 2.4, the I p evolution is sensitive to the temperature evolution, determined by the parameters T initial , T CQ and t TQ . T initial can be measured and is thus taken from experimental data, but t TQ and T CQ are regarded as free parameters. The hot-tail RE generation is exponentially sensitive to t TQ [28], with a quick TQ (small t TQ ) resulting in a large hot-tail RE generation early in the disruption and consequently a higher post-CQ I p . An increased t TQ results in a smaller RE generation, see figure  2(a). The hot-tail seed loss parameter can obviously be adjusted to counteract the effect of a small t TQ , but a small t TQ results not only in increased hot-tail RE generation, but also increased Dreicer RE generation, so the effect of a small t TQ can not be completely cancelled by assuming that the hot-tail seed is lost. As shown in figure 2(b), variations of f HT below 1% have a minor impact on the results for the typical size of the hot-tail seed in these simulations, because at this point, Dreicer generation becomes the dominant RE generation mechanism.
The parameter T CQ affects the I p decay rate during the CQ. Since the duration of the TQ (∼1 ms) is much shorter than the duration of the CQ (∼5 ms), T e ≈ T CQ for most of the CQ and the plasma resistivity σ ∝ T 3/2 e determines the I p decay rate, see figure 2(c). We can therefore infer T CQ by choosing it to match the initial I p decay rate. Note that T CQ is not necessarily the prevailing temperature during the RE plateau phase, but only during most of the CQ, i.e. the temperature may well decay further later in the disruption. The black dashed line in figure 2(c) shows the current evolution assuming a flat temperature profile with a value corresponding to a homogeneous spreading of the initial thermal energy (T e = 1.24 keV). Such an assumption leads to 30% higher final runaway current, since the runaway generation becomes more efficient in a larger volume, but the overall evolution of the current remains qualitatively the same.
The amount of Ar assimilated in the plasma, quantified by the parameter r Ar/D , has a direct effect on I p , a higher r Ar/D leading to a lower post-CQ I p , see figure 2(d), until the post-CQ I p approaches zero.
The search of the parameter space was done iteratively, starting from the set of parameters used in the initial hybrid temperature model simulations. The values of r Ar/D used in the initial simulations were kept the same in these new simulations. T CQ was then chosen to match the I p decay rate during the CQ, defined as ∆I p /∆t during half of the CQ (while I p decreased through 25% to 75% of the CQ). Then the hot-tail seed loss was adjusted to match the measured post-CQ I p . For discharges #34183 and #35649, the post-CQ I p remained too high even when all hot-tail RE seed was removed, so t TQ was increased, keeping the hot-tail seed loss at 100%, until the the measured post-CQ I p was matched. Since t TQ also affects the I p decay rate to a small extent, the process had to be iterated a few times. The conditions for matching were that ∆I p /∆t should be within 10% of the measured value, and that the post-CQ I p should be within 10 kA from the measured value. The chosen values of the simulation parameters are listed in table 5.
Again using discharge #34183 as an example, the I p evolution and the RE generation rates are shown in figure 3. The corresponding figures for the other three modelled discharges are similar. The transition between the CQ and the RE plateau phase is less marked in these simulations than when the hybrid temperature model was used, since the TQ ends less abruptly when prescribed by equation (2), as seen when comparing figure 3(a) with figure 1(a). The exponentially decaying temperature does not represent the physical T e evolution exactly, but is a better approximation than the significant re-heating predicted when switching to the time-dependent energy transport equations as was done with the hybrid temperature model. The RE generation rates shown in figure 3(b) differ from those shown in figure  1(b) mainly in the almost completely suppressed hottail generation.
When the hot-tail generation is suppressed, the electric field is allowed to build up, resulting in larger Dreicer generation.

JET simulations with the hybrid temperature model
Simulations of the JET disruptions were first attempted with the hybrid temperature model, but with similar complications as for the AUG discharges, i.e. for discharges #92454 and #92460 no combination of parameters could be found to match the measured I p evolution. For JET discharges #95125 and #95129, the energy transport calculations predicted a continued drop of the on-axis temperature from T switch = 100 eV down to approximately 20 eV where it remained for a ms after continuing down to the final equilibrium temperature of 1 eV. This is shown for #95125 in figure 4, where the simulation parameter values, found through an iterative search similar to that described in section 4.2, were t TQ = 0.175 ms and r Ar/D = 0.25. The resulting I p evolution and RE generation rates are shown in figure 5, and, as shown, the I p evolution was reproduced with this temperature evolution. It was concluded that, again, the failure to reproduce the I p evolution in the other discharges was mainly a consequence of a failure to predict the temperature evolution. Hence, simulations were performed with a prescribed exponential temperature drop down to a temperature of approximately 20 eV, which made it possible to reproduce the I p evolution also for the other JET discharges listed in table 3.  Figure 5: (a) Current evolution and (b) RE generation rates for JET discharge #95125 using the hybrid temperature model with a switch to energy transport calculations at 100 eV. Note that after ∼ 29 ms, the RE generation is dominated by the avalanche mechanism.

JET simulations with the exponential temperature model
As stated in section 3.2.2, the fraction of the injected Ar atoms that assimilate in the plasma, f Ar , is constrained by the condition that the maximum calculated line integrated free electron density, which is very sensitive to the assimilated Ar density, should not deviate by more than 10% from the measured value. Since the calculated free electron density is not very sensitive to the other simulation parameters, the Ar assimilation fraction, and thereby the ratio r Ar/D , was chosen first.
The other parameters were found iteratively, as for the AUG discharges.
Once the assimilation fraction f Ar , and hence the values of r Ar/D (r Ar/D = f Ar N Ar /(V p n e0 ) = 0.5N Ar /(V p n e0 )), were chosen, T CQ was chosen to match the I p decay rate, and the value of t TQ was chosen to match the measured post-CQ I p . In all cases, since Dreicer was the dominant RE generation mechanism, the I p evolution could be modelled without assuming hot-tail losses, so the hot-tail seed loss was neglected. The conditions for matching were the same as for the AUG discharges, i.e. that ∆I p /∆t (for definition see section 4.2) should be within 10% of the measured value, and that the post-CQ I p should be within 10 kA from the measured value.
The I p evolution shown in figure 6 shows a less abrupt transition from the CQ to the RE plateau phase than the measured data, just as for the AUG case in figure 3(a). This is, once again, due to the smooth end of the TQ given by the exponential approximation. The comparatively small hot-tail peak is visible to the left in figure 6(b), whereas the Dreicer generation continues for several ms. This differs significantly from the generation rates shown in figure 5(b), where the Dreicer peak is almost as narrow as the hot-tail peak.
The calculated line-integrated free electron density is shown together with a measurement-based estimate of the same quantity in figure 6(c). The estimate is based on the line integrated density measured by interferometry (KG1 diagnostic, corrected for fringe jumps). The much faster increase of the calculated density immediately at the start of the simulation is due to the assumption that all Ar enters the plasma instantaneously, which is obviously not correct, but has a minor impact on the calculated plasma current.
Using only the exponential temperature model, the free electron temperature profile remains the same throughout the simulation, as shown in figure 7(a). The free electron density profile is re-calculated to be consistent with the temperature profile, as the ionization states of the injected argon changes with the temperature. Since the radial distribution of the injected argon is approximated by the initial distribution of the deuterium, the density profile remains nearly constant throughout the simulation, as shown in figure 7(b). Some small changes in the profile shape can be seen when comparing the density profile at 26.9 ms and at the end of the simulation (60.0 ms).

Discussion and conclusions
The most important parameters for the simulated discharges in JET and AUG are listed in table 1 and the simulation parameters in table 5. We note that the initial free electron densities and also the externally applied magnetic field are similar for the two sets of discharges, whereas the machine size and the initial plasma current are significantly larger in JET.
It should be noted that the higher current in JET is partly a consequence of the larger machine size -the current densities are similar (1.5 -1.8 MA/m 2 for AUG and 1.9 -2.0 MA/m 2 for JET). In AUG, the initial temperature and the injected Ar density (r Ar/D · n e0 ) are larger than in the JET discharges. The simulation parameters needed to reproduce the I p evolution when using the same modelling strategy are similar, except that the hot-tail seed survival fraction had to be chosen very small to reproduce the AUG discharges, whereas its value was unimportant for modelling of the JET discharges, due to the rather small hot-tail generation. The argon assimilation rates (expressed through the argonto-deuterium density ratio r Ar/D ) are based on experimental data as discussed in section 2.3, and are consistently larger in the AUG discharges.
Our simulations indicate that the hot-tail RE generation is much smaller in JET as compared with AUG. This may be understood by the fact that at a fixed TQ time, the hot-tail seed increases exponentially with initial temperature due to the longer slowingdown time of the hot-tail electrons. According to the simulations, the TQ times in the two devices appear to be similar, although this property would be expected to scale with machine size [1]. This may be because the higher initial temperature in the AUG discharges is compensated by a relatively higher density of impurities (due to the smaller plasma volume), which can then radiate the heat more efficiently.
An important remaining problem to be solved before reliable predictive simulations can be made is the self-consistent modelling of the temperature evolution. Although we found that a time-dependent energy transport model, including ohmic heating and impurity radiation losses, allowed for reproduction of the plasma current evolution for some JET discharges, this model was prone to predict a re-heating of the plasma after a forced temperature drop in other discharges in a way that contradicts the experimental data.
This suggests that some additional losses are present, such as radiation from wall impurities and transport losses remaining during the CQ. The importance of transport losses is expected to be higher in a smaller machine, which might explain why the prediction of an experimentally excluded reheating was less common for JET than for AUG. Time dependent radial profiles of the magnetic field fluctuations are difficult to obtain from experiment, although their amplitude during the current quench has often been found to be fairly small [43]. Numerical modelling of the transport due to magnetic fluctuations [44] would therefore require exploration of a fairly large additional parameter space and the effect of this on the results presented remains an open question.
The parameter T CQ used in the exponential temperature model lies consistently around 20 eV for (a) (b) (c) Figure 6: (a) Current evolution, (b) RE generation rates and (c) line-integrated free electron densities for JET discharge #95125 with exponential temperature decay to 21 eV. The reason for the decrease in the measured density after ∼40 ms is probably that the actual temperature has dropped to the ∼1 eV range in the absence of an ohmic current and associated ohmic heating, leading to lower ionization. both machines. It is assumed that whereas this temperature prevails during most of the CQ, it is only maintained during this period through collisional heating by the remaining ohmic current, and the temperature during the RE plateau is likely to be about 1 eV on average. The existence of such a temperature plateau is predicted by the energy transport equations in the cases when the I p evolution is well reproduced by this model. Experimental data can neither confirm nor exclude the existence of such a temperature plateau.
The main conclusion of the presented set of simulations is that the RE generation models in go are able to reproduce experimentally measured runaway currents for both AUG and JET using similar assumptions and modelling strategies. This increases our confidence that the RE models adequately describe the runaway physics in tokamaks of different sizes. The relative simplicity of the code makes it possible to run a large number of simulations and qualitatively investigate the consequences of the simultaneous variation of several different parameters over large intervals.