Time-dependent runaway electron simulations: Ampere–Faraday equations implemented in CQL3D

The runaway electron (RE) distributions driven by a large toroidal electric field induced by the drop in the temperature profile due to disruption or pellets are comprehensively simulated by the 3D Fokker–Planck (FP) solver CQL3D (Harvey and McCoy 1992 Proc. of IAEA TCM), recently coupled to the Ampere–Faraday (AF) equations. The evolution of the toroidal current in a plasma occurs on a resistive time scale, τres  =  4πa2/(c2η), which is typically of the order of seconds in present tokamaks. Here, a and η are respectively plasma radius or radial extent of a current density perturbation, and Ohmic resistivity. From the Faraday EM equation, the toroidal electric field is proportional to the time derivative of the poloidal magnetic field, which, from the Ampere equation, is proportional to the toroidal current. Thus, the toroidal electric field rapidly increases due to an abrupt temperature drop decrease in conductivity, to prevent change in the toroidal current faster than τres. This is a example of Lenz’s law. For example, in simulations with KPRAD (Whyte et al 2003 J. Nucl. Mater. 313–6 1239) of neon pellet injection into a DIII-D shot, Te drops from 2 keV to 10 eV in 0.1 ms and Zeff increases 1–4, giving that the electric field increases 3500×  to 0.8 V cm−1. As described in Harvey et al (2000 PoP 7 4590), this places much of the tail electron distribution beyond the Dreicer runaway velocity, giving so-called ‘hot-tail runaways’ which for a time are the dominant source of runaways, more so than the knockon source. In this prior calculation, performed for a single flux surface, the toroidal current density is held constant, on the basis that τres is large. Most of the initial current can be converted to runaway current, which is then dangerous, particularly for ITER. A more comprehensive A–F model recently implemented in CQL3D, taking into account the time-development of the full-plasma-width toroidal electric field on time-scales of order τres applies an iterative technique for the toroidal field previously developed for a different application (Kupfer et al 1996 PoP 3 3644), maintaining the implicit-in-time evolution of CQL3D. The degree of runaway current formation is reduced in AF augmented CQL3D, but the basic mechanism of ‘hot-tail runaways’ remains a dominant contribution to the REs at early times after the Te drop in these simulations. On the other hand, a NIMROD (Sovinec et al 2004 J. Comput. Phys. 195 355) simulation of shattered-pellet shutdown of DIII-D plasma (Kim 2018 APS/DPP Meeting), gives a slower thermal quench; when the plasma profiles and electric field are coupled one-way to CQL3D, the ‘hot-tail’ REs are much less, and growth of RE is dominated by the knockon process.

The runaway electron (RE) distributions driven by a large toroidal electric field induced by the drop in the temperature profile due to disruption or pellets are comprehensively simulated by the 3D Fokker-Planck (FP) solver CQL3D (Harvey and McCoy 1992 Proc. of IAEA TCM), recently coupled to the Ampere-Faraday (AF) equations. The evolution of the toroidal current in a plasma occurs on a resistive time scale, τ res = 4πa 2 /(c 2 η), which is typically of the order of seconds in present tokamaks. Here, a and η are respectively plasma radius or radial extent of a current density perturbation, and Ohmic resistivity. From the Faraday EM equation, the toroidal electric field is proportional to the time derivative of the poloidal magnetic field, which, from the Ampere equation, is proportional to the toroidal current. Thus, the toroidal electric field rapidly increases due to an abrupt temperature drop decrease in conductivity, to prevent change in the toroidal current faster than τ res . This is a example of Lenz's law. For example, in simulations with KPRAD (Whyte et al 2003 J. Nucl. Mater. 313-6 1239) of neon pellet injection into a DIII-D shot, T e drops from 2 keV to 10 eV in 0.1 ms and Z eff increases 1-4, giving that the electric field increases 3500× to 0.8 V cm −1 . As described in Harvey et al (2000 PoP 7 4590), this places much of the tail electron distribution beyond the Dreicer runaway velocity, giving so-called 'hot-tail runaways' which for a time are the dominant source of runaways, more so than the knockon source. In this prior calculation, performed for a single flux surface, the toroidal current density is held constant, on the basis that τ res is large. Most of the initial current can be converted to runaway current, which is then dangerous, particularly for ITER. A more comprehensive A-F model recently implemented in CQL3D, taking into account the time-development of the full-plasma-width toroidal electric field on time-scales of order τ res applies an iterative technique for the toroidal field previously developed for a different application (Kupfer et al 1996 PoP 3 3644), maintaining the implicit-in-time evolution of CQL3D. The degree of runaway current formation is reduced in AF augmented CQL3D, but the basic mechanism of 'hot-tail runaways' remains a dominant contribution to the REs at early times after the T e drop in these simulations. On the other hand, a NIMROD (Sovinec et al 2004 J. Comput. Phys. 195 355) simulation of shattered-pellet shutdown of DIII-D plasma  APS/DPP Meeting), gives a slower thermal quench; Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Runaway electron (RE) distributions in tokamaks can be induced by the large toroidal electric field resulting from a sudden drop in the temperature profile due to disruption, pellet injection, or massive gas injection. Such REs can be extremely dangerous to the machine vacuum vessel, and are considered a major, only partially resolved issue for the ITER tokamak. The CQL3D Fokker-Planck (FP) code provides a comprehensive computational model for time-dependent bounce-averaged (BA) electron and ion distribution functions in tokamaks, and in particular for REs and their interactions with radiowaves [1]. In this paper, we present first results of coupling the FP code to the Ampere-Faraday (AF) equations for the self-consistent, time-dependent toroidal electric field.
In the present paper, application to RE modeling is made using the zero-orbit-width (ZOW) version of the finite-difference continuum CQL3D. Averaging the relativistic Boltzmann kinetic equation over gyro-angle, toroidal angle, and bounceperiod gives the 3D evolution equation for plasma distribution function to be solved: where f 0 is the electron distribution as a function of (u 0 , θ 0 ), the momentum-per-mass and pitch angle at the minimum-B point on each flux surface, ρ the generalized radial coordinate, and t the time. Jacobian λ = v ||0 τ b accounts for the volume of the flux surface occupied per perpendicular area at the midplane. v ||0 is parallel speed, and τ b the bounce period. The midplane quantities u 0 , θ 0 , ρ form the constants of motion for the ZOW orbits given by u(s) = u 0 , sin 2 (θ(s)) = (B(s)/B(s = 0))sin 2 (θ 0 ), where s is distance along a magnetic field line measured from the minimum-B point. The Liouville theorem gives that f (u(s), θ(s), ρ) = f 0 (u 0 , θ 0 , ρ). The . . . represents bounce-averaging (BA) along the trapped or passing orbit in the tokamak given by the COM equations [1,2]. Operators R and S represent a divergence free radial diffusion and pinch operator, and particle sources including the large-angle scattering knockon operator, and velocity-dependent sinks. Operator C is the fully relativistic BA nonlinear collision operator [3], although usually the quasi-relativistic operator is sufficient [1,4,5]; Q is the full BA Stix QL relativistic operator [6]; E is the BA toroidal electric field term; and H is the BA synchrotron and bremsstrahlung radiation term. Units are cgs Gaussian, except convenient units are used in the figures.
These equations were used in [7,8] to demonstrate on a single flux surface that 'prompt runaway' could dominate the production of REs during the thermal quench phase, and for some time afterwards, as opposed to Dreicer 'drizzle' and knockon produced REs. The thermal quench was based on KPRAD modeling [9] of n e , T e , T i , Z eff time-dependent profiles. Also, relatively small levels of stochastic magnetic field induced radial diffusion could terminate buildup of large RE populations. In carrying out these calculations, the plasma current density was kept constant in accord with Lenz's Law, by iteratively adjusting the toroidal electric field.
The bulk, low-velocity portions of the distribution functions are maintained in sync with the background densities and temperatures from, for example, KPRAD, by a 'partially nonlinear' procedure [1,7,10] applied to the P 0 -expansion term in the Legendre expansion of the f e u 0 -dependence. That is, the collision coefficients are obtained at each radius from a Legendre polynomial expansion of the f e in pitch angle. The partial linearization refers to resetting the P 0 expansion term, the pitch angle average of the distribution at each u 0 , to a specified Maxwellain distribution. This simulates the effects of transport processes on the bulk, low-velocity plasma, but does not change the plasma current carried by the distribution.
In the present paper, the previous calculations [7] are extended to full plasma radius. With the new AF addition to CQL3D, the current is allowed to relax self-consistently by solving the FP equation simultaneously with the toroidal electric field for the AF equations. This enables comparison with the prior single flux surface runs. For comparison purposes, a background plasma evolution is assumed, applying the single surface time-dependence [7] at all plasma radii. As a result, the toroidal current is found to be reduced considerably during and immediately after the thermal quench. Nevertheless, the 'prompt runaway' process remains a dominant source of RE seed current, and leads to a current quench phase. Small levels of magnetic turbulence are required in order to terminate the RE buildup.
In the following, the AF equations are given, and the method for obtaining a time-implicit solution simultaneously with the CQL3D code is indicated. To the authors' knowledge, this feature is new in continuum FP modeling. A simultaneous solution for RE kinetics with AF equations was previously investigated with the ARENA code using a Monte Carlo approach for the electron modeling; the authors stressed that the simultaneous solution was crucial for accurate RE simulations [10], that is, the constant current density assumption such as used in [1] may not be well satisfied.. In the following, a verification of the combined CQL3D-AF solution is given, as compared to a simple constant-current simulation of [1]. That is, a time-dependent, full-radius simulation of thermal-and early current-quench phase is given. In addition, preliminary results are given for one-way coupling of the NIMROD MHD code [11,12] to CQL3D, and subsequent knockon induced avalanche.
In summary, two new aspects of RE simulations reported here are: (1) The calculation of a simultaneous solution of the 2 V-and-radius continuum BA FP RE distributions with CQL3D including, crucially, the simultaneous solution of the AF equations for the toroidal electric field; and (2) first coupling of CQL3D to the NIMROD time-dependent MHD code, providing an alternative calculation of the toroidal electric field.

Solving AF equations for the toroidal electric field
The evolution of the toroidal current in a plasma occurs on a resistive time scale, τ res = 4πa 2 /(c 2 η), which is typically of the order of seconds in present tokamaks. Here, a and η are respectively plasma radius or radial extent of a current density perturbation, and Ohmic resistivity. From the Faraday electromagnetic equation, the toroidal electric field is proportional to the time derivative of the poloidal magnetic field, which, from the Ampere equation is proportional to the toroidal current. Thus, the toroidal electric field increases due to rapid temper ature drops, to prevent change in the toroidal current faster than τ res . For example, in simulations with KPRAD of neon pellet injection into a DIII-D shot, T e drops from 2 keV to 10 eV in 0.1 ms and Z eff increases from 1 to 4 [7], giving that the electric field increases 3500× to 0.8 V cm −1 . As described in [7], this rapid cooling of the bulk electrons results in much of the less-collisional tail electron distribution being beyond the runaway velocity, giving so-called prompt 'hot-tail runaways' which for a time are the dominant source of runaways, more so than the knockon source and Dreicer 'drizzle'. In this prior calculation, performed for a single flux surface, the toroidal current density is held constant, on the basis that   τ res is large. Most of the initial current can be converted to runaway cur rent. This process is described more comprehensively by the AF equations coupled to CQL3D. The Faraday and Ampere equations written in toroidally symmetric geometry can be cast, respectively, as Where V loop is loop voltage at each flux surface obtained from the edge value and time-derivative of the poloidal flux between the plasma edge and ρ, and the poloidal magnetic field is given by the integral of the toroidal current up to radius ρ. Using RB φ is constant on a flux surface, and that the flux surface shapes are constant in time, E φ (ρ,θ pol ) = E 0 (ρ)R mag /R, then equations (3) and (4) can be written in terms of B θ0 , J ||0 evaluated at the outer midplane of the flux surfaces. The quanti ty J ||0 is obtained from the parallel velocity moment of the CQL3D distribution f 0 . A simultaneous fully implicit solve of the combined equations is computationally very difficult, since the 2D momenta and 1D radius are all coupled through the AF equations, giving typically a 10 6 square matrix to be inverted which is very large. On the other hand, the explicitin-time solution for E 0 is numerically unstable [11]. An iterative method has been used, suggested by an approach developed by Kupfer [14] in a somewhat different context. The advanced time electric field is expanded as E n+1 = E n + ΔE where in general ΔE will be small, and similarly the time advanced distribution function is expanded proportional to ΔE. Amazingly, this divides solutions of the FP equation and the AF equations into separate time-implicit solves, the latter being tridiagonal. Typically, two iterations are sufficient per time-step for convergence. This gives only a factor ~2 increase in the computer time required for time stepping the combined 3D FP-AF equations.
A verification test has been constructed for the combined CQL3D-AF comparing the computed total current and onaxis electric field from the space-time solution with analytic solutions for normal modes of AF in a uniform cylinder. For uniform cylindrical geometry, the AF equation reduces simply to For BC E z (r = a, t) = 0, the λ k are the kth zeros of the Bessel function. Thus solutions J 0 (λ k r) will exponentially decay with time constant τ = 4πσ/(c 2 λ 2 k ). The straight cylinder is approximated in CQL3D by a large, but not infinite aspect ratio torus (R/a = 1400). Conductivity σ is the Spitzer value, and in CQL3D calculated from solution of the FP equation. In figure 1, a comparison between the analytic and CQL3D solution is shown for decay of the k = 3 mode (i.e. with A 3 = 1 and A j = 0 for j not equal to 3). Figure 1(a) shows the space-time variation of E z from CQL3D, with an insert showing convergence of the solution in two steps. Figure 1(b) gives the CQL3D-AF time-decay of the central E z (r = 0,t) (the dotted solid line) in comparison with the analytic decay expression using conductivity according to Spitzer (upper, green curve), and according to CQL3D at r = 0 and at r = a. There is only a small variation of conductivity versus r/a. The central CQL3D conductivity agrees reasonably well with Spitzer, and CQL3D-AF agrees well using its internally calculated conductivity. This is taken as a good verification of the CQL3D-AF code.

RE calculation due to pellet-induced temper ature drop across the whole plasma radius: AF E(ρ, t) solution
In [7], prompt 'hot-tail' REs remaining when the plasma cools faster than the tail electron collision time, were identified as a major source of REs, either as the seed for a knockon avalanche, or the main source of REs. The electric field driving this phenomena was determined by holding the total current density constant during the T e -drop simulation, in general accord with Lenz's Law, an approximation which will overestimate the strength of the induced toroidal electric field. Rather, the current density will decrease for a resistive time, which is becoming shorter as T e decreases, until the highly conductive REs appear. This phenomena was examined on a single representative plasma flux surface. It was also shown that small levels of radial diffusion could cause sufficient radial transport of the REs to dissipate them to the chamber wall or divertor.
The single flux-surface simulations [7], were based on background plasma evolution given by the KPRAD [9] pellet code. The electron temperature dropped from 2 keV to 10 eV in 0.1 ms, and then held steady, as shown in figure 2. Also, during this period, the electron density increased from 7 × 10 13 to 2 × 10 14 cm −3 , and Z eff from 1 to 4. To compare with [7], and as an early application of the CQL3D-AF code, this time-dependence of background T e /n e /Z eff is kept the same across the entire plasma radial profile in the present full-radius modeling. For an initial 1.0 ms before the pellet plasma variation begins, a steady toroidal current at a DIII-D-like value is established through application of a constant in radius and time toroidal electric field in the T e = 2.0 keV plasma, with no simultaneous AF equation solve. Large time-steps, 0.1 ms, are used in this phase, well-beyond the collision time, giving essentially an iteration to near steady-state electron distributions under the influence of the initial applied steady electric field; the procedure is numerically stable at such large time steps in view of the fully time-implicit nature of the numerical solution. The AF equation solve, simultaneous with the rest of CQL3D, is turned on at t = 1.0 ms in the simulation, and the plasma parameters are varied throughout as in figure 2. This is a model problem, for comparison with [7], in which background plasma profiles are constant versus radius. The boundary condition for the toroidal electric field at the plasma edge is that it is held constant at the established low value at t = 1 ms. (Future work will include choosing a boundary electric field consistent with the presence of a conducting vacuum vessel, and evolving plasma equilibrium. One approach is by coupling to an MHD code, as discussed below.) The calculated toroidal electric field is shown in figure 3, and the corresponding total current in figure 4. The large E φ spike is induced as a result of the sharp pellet-induced T edecrease across the plasma radius. The total current decreases by ~20%, with a resistive time ~0.5 ms during the T e -drop, and then, as the prompt 'hot-tail runaway' distribution is established, the plasma conductivity greatly increases and the slow current quench phase begins. We emphasize that the solution is for the whole velocity range of the distribution at each radius, from v = 0, v max , where v max is generally taken to be large enough that there is negligible flux off the edge of the A critical quantity for both the RE process and the possibility of a deleterious knockon avalanche of electrons is the ratio of electric field to critical field for electron runaway, E crit = 4πn e e 3 ln(λ)/(m e c 2 ). The critical momentum-per-mass, above which the electric field force exceeds electron Coulomb drag and electrons thus runaway, is u crit /c = 1/(E/E crit − 1) 1/2 ; this quantity is shown in figure 5 which particularly illuminates the hot-tail runaway phenomena. E/E crit has jumped up to 75 with the T e -drop, and consequently the critical velocity has dropped to ~0.1c in less than 0.07 ms, whereas 4v Te electrons in the initially T e = 2 keV electron distribution have collisional slowing down time ~0.08 ms, giving the 'hot-tail' REs. If more current were lost during the fast T e -drop than in figure 4, then u crit would be larger, and the 'hot-tail runaways' exponentially less and hence would form a smaller runaway population and seed population for a possible knockon avalanche. The formation of a substantial hot-tail runaway population is evidently a very sensitive function of the thermal quench dynamics, and warrants detailed studies. Figure 6 shows cuts through the electron distribution at constant pitch angle on the top row, at three times as the runaway develops, and shows the specific parallel current versus u/c on the bottom row. At the T e -drop, the current becomes almost entirely carried by the prompt 'hot-tail runaway' electrons, as shown in figures 6(a) and (d). The near vertical line at u/c ~ 0.0 is the thermal distribution and thermal current contribution. The clump of REs advects to higher velocity and energy with time, as shown in the remaining subplots. The electric field at the end of the run is ~0.01 V cm −1 . Electrons at speed c, gain energy 0.3 MeV ms −1 , giving du/c = dE/(mc 2 ) = 0.6/ ms, consistent with the acceleration shown in the figure. In the present simulation, the maximum energy, 100 MeV, corresponds to u/c ~ 200, and thus virtually all electrons remain on the grid. These simulation results also reinforce the point first made in [11] with the ARENA code, that a comprehensive simulation of disruption RE production requires the simultaneous solution of the electron kinetic and AF equations.
The space-time evolution of runaway current j run , defined as that portion of the current from electrons at u ⩾ u crit , is given in figure 7. At the larger times shown, j run is eroding at the plasma edge, consistent with resistive time τ res = 4πL 2 / (c 2 η || ) = 0.5 ms, using scale L = 10 cms; and is concentrating slowly towards the plasma center, consistent with τ res = 200 ms.
The robustness of the above RE production against stochastic magnetic field fluctuation radial diffusion has been examined for two values of D rr0 = 1 and 10 m 2 s −1 , applying electron radial diffusion D rr = D rr0 × (v || /v th0 ) across the discharge, in accord with Rechester and Rosenbluth [15]. These values correspond to magnetic fluctuation levels δB/B = 0.04% and 0.08%. Both cases result in dissipation of the RE discharge, in 2 ms and 9 ms respectively. That is, it still holds [7] that small levels of stochastic fields, if they extend throughout the radial cross-section, can empty the REs to the plasma edge. However, some modeling indicates only limited radial extent of the stochastic regions, and only for a limited post-thermal-quench period [16,17]. This will be addressed in future work, through the coupling with NIMROD.
The knockon source term [7,8] is included in these calculations. From figure 8(a) showing the velocity-space source at radius r = 0.5a at 1 ms after the pellet, about half the knockon source electrons are born as trapped particles. The pitch angle average source [7] in figure 8(b) exhibits the classical exponential fall-off with energy. At the minimum u crit /c = 0.1 just after the pellet, this velocity and the prevailing plasma parameters give pitch angle scattering time τ ei = 1.3 µs which is much less than the RE evolution time. Thus, the pitch angle structure of the knockon source will be smoothed out by collisions. Similarly, at later times with typical u crit /c = 0.4 the scattering time is still small.
In the present simulation, the knockon source of particles is small compared to the prompt 'hot-tail' electrons. The knockon electron growth rate is γ 0 ~ 2(E/E crit ) s −1 in the postpellet plasma, giving at this time in the simulation an inconsequential maximum rate of 140 s −1 ; it was larger, 460 s −1 , in [7], but still inconsequential. Apparently the knockon source will become important in cases of substantially slower thermal quench, such as in the following preliminary CQL3D modeling with NIMROD plasma parameters and electric field. In addition in the present simulation, bremmstrahlung, gyro-, and toroidal-geometry synchrotron radiation are all small.

Computations discussion
In addition to the new, effectively time-implicit, simultaneous solution of the AF equations, several features of CQL3D enable the efficient solution for RE dynamics. Variable grid spacing in momentum-per-mass space facilitates the solution on a single mesh for both the very low temperature post-disruption electrons along with the energetic RE tail to 100 MeV and beyond. A geometric mesh in u is used for this. Geometric packing of the grid in pitch angle is also used, in this case to increase resolution in the parallel direction. Use of a 'quasirelativistic' operator [4], which has been benchmarked against the fully relativistic operator for non-relativistic bulk plasma [5], greatly speeds and simplifies the computations. Computer cpu time per time-step/per radial point, is typically ~1 s. Variable time-step facilitates can be used to keep the time-step shorter than the solution time-scales of interest; alternatively, large time steps provide a means to iterate to steady-state solutions in 5--0 steps. The open-software SPARSKIT [18] sparse matrix GMRES solver is used for the fully-implicit 3D solves requiring only 5.0 h single processor cpu time for 4.5 ms simulation with 350 u-, 80 θ-, and 51 radial-mesh points, and 161 time steps. MPI parallelization over flux surface grid points is implemented for cases without radial diffusion. A parallelized GMRES is being investigated for cases with radial diffusion.
The code is open-source [19].
Because of the crucial nature of controlling REs for the safety of ITER, there are many recent calculations on this topic, for example, [20][21][22]. The major differences between the calculations reported here and in these references are that a 2D-in-v is used for the complete BA electron velocity distribution here, as compared to a lumped density and a 1D-in-v FP distribution of the REs in [20,21] and a 2D-in-v, no-BA, high-velocity limit FP solution in [22]. And, in [20,21], the solution of AF equations for E φ (ρ,t) is accomplished using the more limited lumped-RE/1-v FP solution, whereas in [22], a constant current density approximation, as in [7], is made. The treatment of the thermal quench distribution dynamics reported here is thus more complete. On the other hand, [20,21] are coupled to equations for the radial transport of plasma energy, whereas the calculations reported above are based on the 0-radial-D KPRAD solution extended across the plasma radius. Aleynikov and Breizman [22] includes important scattering of REs on screened high-Z nuclei [23][24][25], not  Coupling of CQL3D with a NIMROD shattered pellet simulation gives an exponential growth of RE density and current in general accord with the knockon growth rate. (E/E crit ) ~ 500 during the growth phase, giving growth rate ~1/s. The green curves are at r = 0.75a, the blue at 0.52a, and the red at 0.25a, and the start-times reflect the progress of the pellet into the plasma. Prompt 'hot-tail' REs provide an initial knockon current growth, although is this case the initial RE current density is a very small fraction of the total. yet added to CQL3D, but particularly important in the current quench phase. By comparison, CQL3D includes important radial transport of the whole distribution function. For these several RE models, there is also a need for validated models for coupling to space-time dependent impurity deposition from disruption or mitigation scenarios.

One-way coupling of nimrod to CQL3D
The CQL3D-AF model has not (yet) included motion of the equilibrium. This is a crucial feature for RE modeling which is being addressed by coupling to the NIMROD 3D extended MHD and time-dependent transport code [12], including a shattered pellet deposition model [13]. Thus far, an initial one-way coupling has been accomplished, that is, NIMROD provides plasma profiles and toroidal electric field [12,13], and these quantities have been used to drive the CQL3D code. The NIMROD code has included the AF equations assuming Spitzer plasma conductivity, but does not include REs which increase the conductivity and would slow down the thermal quench reduction in plasma current and give an extended current quench. In the detailed Kim [13] fractured Ne pellet modeling, the thermal quench time τ quench is reduced from the KPRAD derived time ~0.1 ms shown as in figure 2, to 2.5 ms. As shown in figure 9, from CQL3D with NIMROD profile data, this gives very small prompt 'hot-tail' electrons, as expected, given the sensitivity of this process. The E/E crit ratio is ~500 across most of the plasma, giving knockon growth rate ~0.5/ ms at the prevailing density and Z eff . The central runaway current is 10 −3 of the total current density, 10 −5 at r = 0.5a, and 10 −8 at r = 0.75a, indicating that the knockon avalanche will concentrate the current in the plasma center. When the CQL3D runaway current is coupled back to NIMROD, it will slow down the plasma current evolution, giving the current quench.

Summary and conclusions
This paper has presented a new, comprehensive continuum BA nonlinear 3D FP computation of RE dynamics including crucial CQL3D coupling with the AF equation for the self-consistent toroidal electric field, and radial diffusion and pinch of the distribution. Both the thermal and slower current quench RE phases are modeled. The code is computationally efficient, and affords a basis for broad investigation of RE scenarios. A particular sensitivity of the fraction of plasma current drop during the thermal quench, to the rate of the thermal quench is indicated. First results of a one-way coupling of time-dependent data from the NIMROD MHD code into CQL3D are also given, and show much slower thermal quench than was obtained with the original KPRAD based modeling. Improved models of the physics of the thermal quench (disruption, pellet ablation, massive gas injection) will be required to accurately estimate the resistive drop in the plasma current during the thermal quench; the heat energy is thereby radiated by plasma impurities rather than the magnetic field energy going into the REs, and optimization of this effect will mitigate RE machine damage.