Disruption mitigation modeling for the SPARC tokamak Open Access

The high plasma pressure, large plasma current, and deuterium-tritium operation in the SPARC tokamak¹ will create a unique proving-ground for disruption mitigation strategies that could be deployed in future reactors. In the event that confinement is suddenly lost, the high-field, compact design poses challenges for structural forces and thermal energy mitigation. Similarly, the combination of nuclear seed sources (tritium decay, Compton scattering)² and high avalanche gain (⁠ $∼109$⁠)³ in a high current, D-T tokamak creates runaway electron (RE) mitigation challenges of an order yet unrealized in existing tokamaks. Demonstration of disruption-resilient operation on SPARC could, therefore, be a watershed in the disruption challenge for tokamak reactors.

Two primary tools in the SPARC disruption mitigation strategy are massive gas injection (MGI) and a runaway electron mitigation coil (REMC).⁴ MGI, with one or two injection locations, has been tested on a large number of tokamaks, including DIII-D,^5–8 Alcator C-Mod,^9,10 KSTAR,¹¹ JET,^12,13 ASDEX-Upgrade,¹⁴ and others. More recently, disruption mitigation experimental and modeling efforts on numerous tokamaks have focused on shattered pellet injection (SPI),^15–22 in large part because this technique was chosen for the ITER disruption mitigation system (DMS).^23–25 Unlike MGI, in which the impurities propagate down the injection tube as a gas, SPI accelerates a cryogenic pellet down the injection tube and shatters it into a mixture of fragments, liquid, and gas using a bend or shatter-plate just before reaching the plasma edge. The possible deeper penetration of solid fragments into the edge is one potential advantage of SPI, but the required length of the injection tube on ITER engenders the most significant drawback of MGI. Without a large diameter tube and unacceptably large gas plenum, the gas dynamics would result in extremely slow gas delivery to the plasma over such a distance.²⁵ On the more compact SPARC device, no such obstacle to MGI exists, and the simplicity of MGI (relative to SPI) offers advantages for robustness.

Modeling of MGI with 3D MHD codes has also been carried out for several tokamaks, beginning with NIMROD modeling of Alcator C-Mod²⁶ and then DIII-D.²⁷ Modeling of MGI in JET was also performed with the JOREK code,^28,29 and MGI modeling for SPARC was previously performed with M3D-C1.³⁰ Some notable physics effects seen in MGI simulations include: the more rapid parallel spreading of impurities on rational flux surfaces, due to faster thermal equilibration;³¹ the preferential parallel spreading of impurities toward the high-field-side, due to the “magnetic nozzle” effect;³¹ and the dominant effect of the m = 1/n = 1 convective heat flux on the radiation asymmetry during the thermal quench (TQ).^31,32

In contrast to MGI, the REMC concept^33,34 is untested on a disruption-runaway producing tokamak (although a recent REMC installation on the small HBT-EP device has enabled some physics and engineering validation studies³⁵). The influence of 3D fields on disruption REs has been studied on several tokamaks,^36–39 but there is no well-established solution that is predicted to suppress the runaway electron avalanche in a 9MA, D-T tokamak.⁴⁰ The n = 1 REMC that was designed for SPARC⁴ has been modeled with the NIMROD code,^41–44 and in combination with the ASCOT5⁴⁵ and DREAM⁴⁶ codes, it has predicted full suppression of an otherwise $>5MA$ RE beam.⁴¹ However, modeling also predicts that the n = 1 coil maintains resonance with the plasma only so as long as the safety factor (q) on the axis remains below two, and that good flux surfaces promptly re-form once $q0>2$⁠.⁴⁴ At what point in the CQ this happens is dependent on multiple factors, including the TQ scenario and the simulation boundary conditions.⁴⁴ 

In the following sections, we will first present modeling of MGI in SPARC, in several configurations, including the particular 6-valve layout that was selected for the final design. These simulations predict very high radiated energy fractions and fairly low toroidal radiation peaking for the planned system. In Sec. II, we will present the first modeling that combines the MGI and the REMC applied fields to enable a comprehensive simulation of a SPARC mitigated disruption. The effects of a “closed-circuit” coil on the TQ dynamics are also explored, and the coil perturbation is found to improve the radiated power fraction while also slightly increasing radiation peaking. These simulations again show the importance of maintaining $q0<2$ for the coil operation and highlight some of the sensitivities of the q-profile evolution to simulation parameters. The results point to model improvements that could enable better predictions, which is discussed in Sec. IV.

NIMROD⁴⁷ is a 3D nonlinear extended-MHD code that employs a finite element representation of the poloidal plane and a finite Fourier series in the toroidal direction. The impurity model that is employed in disruption mitigation simulations is described in detail elsewhere.³² In brief, ionization, recombination, and radiation of an impurity species (here neon) are calculated, and the resistive MHD equations are modified to include: a separate continuity equation for each impurity charge state (using a single center-of-mass fluid velocity); energy sinks associated with impurity radiation and ionization; and a local $Zeff$-dependent plasma resistivity. Impurities are deposited as neutrals that have zero temperature; neutrals, both impurities and deuterium, do nothing except diffuse and ionize and have their own diffusion coefficient. Neutral density (of any species) inside the separatrix tends to be very low. Once ionized, they equilibrate (instantly) with the background plasma with a corresponding dilution cooling of the plasma and begin to move with the fluid flow.

In the simulations that follow, all injected neutrals are deposited in the region between the simulation boundary and the equilibrium separatrix. Radially, the deposition profile is flat across this region, with a sharp tanh transition at the separatrix. The deposition patterns have a Gaussian shape in the poloidal direction and a von Mises (periodic normal) distribution in the toroidal direction, with all simulations having multiple injection locations (4 or 6) representing multiple gas valves. An example of a 6-valve deposition profile is shown in Fig. 1, where three upper and three lower valves are modeled with the two rows evenly spaced toroidally and aligned with each other.

The simulations also require the enforcement of a floor temperature (to avoid the occurrence of negative temperatures), which is set to 4 eV and is also the initial temperature in the region outside the separatrix. At any temperature above 4 eV, ionization and radiation will occur at the rate prescribed by the local temperature, but at the floor temperature, no energy is available for these processes since the plasma is not permitted to cool further. The region outside the separatrix is, therefore, effectively clamped at 4 eV, and the deposited neutrals do not ionize until they diffuse across the separatrix. This floor temperature is well below the peak radiation temperature for Ne⁴⁸ (see also  the Appendix), but recombination would be expected to become significant at lower temperatures, so that this approximation could have some effect on the overall dynamics.

Four possible configurations for the SPARC MGI design were simulated with NIMROD, two having six valves and two having four (see Fig. 2). The 6-valve configurations each comprise three valves above the midplane, evenly spaced toroidally, and likewise 3 below the midplane. In the configuration dubbed 6-valve resonant, the upper and lower sets of valves are aligned. In the other case, called the 6-valve non-resonant, the two rows are clocked $60°$ from each other toroidally. Here, resonance refers to injection by upper and lower valves into the same flux tube. The two 4-valve configurations have the pair of upper valves and the pair of lower valves, each $180°$ opposed. The 4-valve resonant option has the lower valves clocked $120°$ from the upper, and the non-resonant option has them clocked $20°$⁠. For each of these layouts, we assume that a 90%/10% ratio of $D2/Ne$ is injected with $1023$ molecules of $D2$ and $1.11×1022$ atoms of Ne distributed evenly among the valves. For the purpose of the simulations, the total gas quantity is assumed to be delivered at a constant rate over a period of 2 ms. Delivering $2×1022$ Ne atoms to the plasma in $≤2$ ms (i.e., more than is modeled here) is a design requirement for the SPARC DMS, and validation of this is presently ongoing, leveraging prototype test stand data as well as 2D computational fluid dynamics simulations. Although simulations begin once the gas reaches the separatrix, the SPARC MGI system has a transit time from the valve to the plasma of 1.3 ms, and the total interval from disruption prediction to gas reaching the plasma may be as fast as 3 ms and no longer than 10 ms.

An initial set of four simulations was performed with toroidal resolution up to n = 10. The simulations have resistive wall boundary conditions applied to the n = 1–6 components, with the resistive wall parameter $vwall(=μ0ηwall/δwall)=40 m/s$⁠, where $ηwall$ is the wall resistivity and $δwall$ is the wall thickness. An ideal wall is applied to the n = 0 component, stabilizing vertical motion, and also to the n = 7–10 components for numerical stability. Two studies^49,50 calculate the relevant wall time for passive vertical stability on SPARC to be on the order of 50 ms, so that little vertical motion is expected during the much shorter TQ phase. The target SPARC equilibrium for these and all simulations presented is shown in Fig. 3 and has a very flat q profile just above one in the center (with $qmin=1.05$⁠). The following simulation parameters are used for all of the simulations in this paper: Spitzer-like resistivity, $η=3.39×10−4Zeff·min{max{Te[eV],4},2000}−1.5Ω·m$⁠; coefficients of isotropic and parallel viscosity, $ν=500 m/s$⁠; parallel and perpendicular thermal conductivity, $κ∥=108m2/s$ and $κ⊥=10m2/s$⁠; hyper resistivity, $η̃=10−6m4/s$ (no other hyper-diffusion is used). Parameters that vary between simulations in different sections, including poloidal and toroidal resolution and density diffusion (for ionized species), are summarized in Table I. Neutral species (Ne and D) have a diffusion of $10m2/s$ for all cases (which is 2–5 times higher than for ionized species). The cutoff temperatures in the evaluation of resistivity (⁠ $4 eV<Te<2000 eV$⁠) do not apply to any other quantity and are not a very significant approximation since most of the dynamics occur within this temperature range (and at higher temperatures reconnection would likely be on faster than resistive time scales).

These simulations have a toroidal source width that is four times broader than what is shown in Fig. 1. Because the source itself is so toroidally broad, these simulations are not likely to give accurate quantitative predictions about metrics such as toroidal radiation peaking but are useful for making general comparisons among the configurations. A second major approximation in this set of simulations concerns the 90% $D2$ in the source. The very large quantity of $D2$ poses numerical challenges for the code, which are overcome by first assuming that 10% of the source $D2$ will be assimilated into the plasma and then by depositing this quantity of D across the entire plasma at the start of the simulation, discarding the remaining $D2$ so that the source is pure Ne. This approximation leaves the full quantity of radiating Ne in the source to be injected over time from the boundary according to the prescribed source distribution, while also allowing the final D content in the plasma to be consistent with what is expected after the injection is complete. The numerical tractability of the simulation is also increased by the lower initial $Te$ and higher density. With an assumed $2×1022$ atoms of D assimilated into the plasma, the target plasma has a $≈3.7$ times higher density and proportionally lower temperature than the expected SPARC flat-top. It should be noted, of course, that localized deuterium in the source would increase the localized cooling of the edge by dilution, which, in turn, affects the contraction of the current profile and the destabilization of MHD modes. As will be described in Sec. II D, this approximation is partially mitigated in the higher resolution simulations but may affect instability onset in all cases.

A comparison of the time-evolution of global quantities for the four simulations is shown in Fig. 4. In all cases, a pre-TQ lasting 2 ms with a radiated power level rising slowly from 5 to 15 GW is followed by a TQ phase with a narrow flash of radiation peaking between 40 and 55 GW. After the TQ, the Ohmic heating quickly rises to match the radiated power in the CQ phase. The thermal energy evolution shows slight trends in which non-resonant configurations have a shorter TQ time than resonant and 6-valves shorter than 4. The peak radiation is unsystematic with respect to these factors. The most obvious variation among the simulations is the different behavior of the 4-valve resonant simulation just after the TQ, when the other three cases exhibit a large plasma current spike. The cause and major consequence of the missing $Ip$-spike for this case is exhibited in the evolution of the n = 1 mode energies and in the radiated energy fractions [Figs. 4(d) and 4(e)]. The large n = 1 mode that peaks just as the current spike begins in three cases does not occur in the 4-valve resonant simulation, where the n = 1 magnetic energy remains more than an order of magnitude smaller at this time. In this case, predominantly even symmetry is preserved until the CQ phase. In addition to the lack of $Ip$-spike, this suppressed n = 1 mode has consequences for the transport of heat to the boundary. At 3 ms, all cases have fully stochastic fields, but the 4-valve resonant simulation maintains field line connection lengths to the boundary that are several times larger than the other cases. In other words, confinement is particularly poor during the $Ip$-spike, allowing fast parallel transport of heat from the core to the boundary, competing with radiation. For this reason, the radiated energy fraction (⁠ $Frad$⁠) is seen to drop to varying degrees as the n = 1 mode peaks, but the 4-valve resonant simulation shows essentially no such drop. These differences in $Frad$ are a matter of a few % however, and all predicted radiation fractions are high.

Because of the toroidally broad Ne injection profiles, the radiation toroidal peaking factor (TPF) is expected to be underpredicted in these cases, especially in the pre-TQ phase. However, a comparison of all four configurations was performed only at this lower resolution, so we present the values of TPF obtained for relative comparison purpose in Table II. Here, we calculate the peaking of radiated energy, integrated over three different phases of the simulation (with times indicated). The TPF is defined as the toroidal maximum over the toroidal average. The differences in the pre-TQ conform to expectations. For instance, the 6-valve non-resonant case has 6 evenly spaced toroidal injection locations and the lowest TFP in that phase. The 4-valve non-resonant case with all injection locations within a $20°$ range has the highest. However, we also find the lowest TPF for the 6-valve non-resonant case in the other phases, with a larger discrepancy than what is seen in the pre-TQ. The 6-valve non-resonant case was selected for the SPARC DMS design, and these low resolution simulations point to one potential advantage. This choice was influenced by the physics understanding at the time of the decision, informed by the simulations reported herein and complementary simulations from the M3D-C1 code,³⁰ together with engineering considerations regarding SPARC port allocations, diagnostic implications, and space constraints. Although a different set of simulations was completed with M3D-C1, it should be noted that for cases with 6-valves, the two codes show very similar levels of peak radiated power in the TQ and comparably high radiated energy fractions.

Only the two 6-valve configurations were modeled at higher resolution. Specifically, the toroidal width of the injection source was reduced by a factor of four to correspond to the distribution shown in Fig. 1. In order to resolve this toroidally localized source, the simulations were run with toroidal mode numbers up to n = 21. The high-resolution simulations were run on twice as many CPUs (1408 vs 704) and took roughly twice as much time to complete (48 vs 24 wall clock hours). Additionally, the approximation that took the deuterium out of the source is partly reversed. That is, in the high-resolution simulations, 20% of the deuterium, for a total of $4×1022$ atoms, is put back into the localized, time-dependent source model, with the remaining 80% of the deuterium assumed to be 10% assimilated already at the start of the simulation. The deuterium in the source is now dominant over the $1.11×1022$ Ne atoms, more closely corresponding to the reality, and the total localized injection quantity increases by almost a factor of 5. The choice of 20% represents an initial effort to restore a substantial fraction of the deuterium to the source; the maximum numerically tractable quantity in the source was not established.

Before looking in greater detail at the high-resolution simulation results, we compare the TQ behavior between the lower and higher resolution cases in Fig. 5. In both cases, we see from the radiated power that with a more localized source (and more deuterium concentrated in the source), the TQ is triggered later, and the radiation peak at the TQ is narrower and larger amplitude, increasing from 50 to 55 GW to 70 GW. The higher resolution simulations also have higher radiated energy fractions than the corresponding low resolution cases. Due to computational expense, both high-resolution cases (magenta) are run just past the end of the TQ until radiation and Ohmic heating balance, but not through the $Ip$-spike phase.

The more localized sources in the high-resolution simulations allow more reasonable predictions for radiation toroidal peaking. The toroidal peaking of radiation, defined as the maximum over the average (along with the toroidal peaking of the Ne quantity), is shown for the two simulations in Fig. 6. Most of the neon remains unassimilated as neutrals outside the separatrix, so the peaking of total Ne mostly reflects the peaking of the source itself, and it can be seen that the non-resonant case with six toroidal locations has about half the TPF for total Ne as the resonant case with three toroidal locations. Although this inevitably leads to lower radiation TPF in the early pre-TQ for the non-resonant case, this does not suggest that the local radiation near each valve is any lower for the non-resonant case. For a more relevant comparison in this phase, the radiation can be integrated over only the upper half of the volume to compare the “upper TPF” for the two cases. For this metric, the two are nearly identical up to 0.25 ms, and the TPF for the resonant case subsequently drops slightly faster. The neon that is assimilated and ionized tends to spread quickly along the field lines and has a toroidal peaking very near to one within 0.5 ms. The peaking of radiation generally exceeds that of ionized Ne throughout the simulation, which can be attributed to various factors at different stages of the simulation, including dominance of lower charge states closer to the source and asymmetric heat flux associated with MHD activity.

The resonant simulation has a clear n = 3 pattern in the radiation in every phase, although at the time of the TQ, the three peaks are very uneven, indicative of a superimposed n = 1 structure associated with the large 1/1 mode. The much flatter radiation pattern at early times in the non-resonant simulations shows a slight n = 6 variation, with the n = 1 being the dominant pattern at the TQ. In both cases, a slight uptick in the radiation peaking occurs at the time of the TQ but not exceeding 1.5. The 1/1 flows that mix neon from the boundary also produce an uptick in the peaking of ionized neon in the non-resonant case.

The relationship between the 1/1 mixing and the TQ radiation flash is evidenced in plots of the radiation pattern and neon density at that time (Fig. 7). Radially, we see that the dominant emissivity at the time of the radiated power peak occurs neither directly near one of the valves nor at the magnetic axis but is maximum at mid-radius and is also poloidally localized. Comparing the emissivity patterns at two toroidal locations (near the toroidal maximum; see Fig. 6), we see a similar structure, rotated poloidally, indicative of a single helical flux tube with the strongest emissivity. The q-profile at this time is very similar to the initial profile with a flat q just above one across much of the core. The plasma has no q = 1 surface at this time, but the radiation peak is at $q≈1.05$ and so is well interior to any low order rational surface such as 4/3. The ionized neon density plots show localized blobs of neon pulled-in from the edge at the same poloidal locations as the peak emissivity. Thus, the peak radiation occurs at the very start of the inward radial mixing of Ne and not when the Ne is well mixed or has reached the axis.

To assess whether SPARC first wall materials can withstand the radiation flash from these simulated high performance mitigated plasmas, the emissivity is ray-traced to the first wall using the Cherab modeling suite⁵¹ (Fig. 8). The peak radiant heat flux exceeds 1 GW/m² but for a duration of only $∼100μ$ s. Calculation of the peak surface temperature using the continuous heat-impact-factor derived in [Hu 2024] finds 20 MJ/m²s^0.5, which is near-to but below the critical value where tungsten heavy alloy (WHA, the material of the main-chamber SPARC limiters) would undergo surface melt and vaporization. It is concluded that SPARC can approach the radiation flash melting limit of WHA in the highest performance pulses, like those modeled here, whereas radiation flash melting of pure W with a critical heat-impact-factor of 50 MJ/m²s^0.5 is unlikely in any SPARC mitigated disruption. It will be an operational goal of SPARC to learn how to reduce this flash melting risk without compromising the high radiated fraction.

One additional analysis was performed for the 6-valve non-resonant simulation, which was to calculate an approximate distribution of radiation as a function of wavelength at the time of the peak radiated power during the TQ. This analysis can be useful in assessing the effect of the radiation flash on materials and is described in the  Appendix.

Initial NIMROD modeling of the SPARC REMC studied only the CQ phase, so that no MHD associated with the TQ was accounted for, and the current profile at the start of the CQ was identical to the flat-top equilibrium profile. That modeling, in conjunction with the ASCOT5 and DREAM codes, predicted full suppression of an otherwise 5MA RE current due to just the REMC perturbations.⁴¹ Subsequent modeling that included a more realistic TQ and associated MHD found that while the TQ itself induced a substantial loss of the RE seed population, the increase in q associated with the current profile rearrangement tended to reduce the resonance of the plasma with the coil in the CQ phase.⁴² The addition of resistive wall boundary conditions to the modeling⁴⁴ further altered the current profile evolution, with the coil proving more effective in the CQ phase in resistive- than in ideal-wall simulations. Building on these results, the REMC modeling is combined here with the MGI model to understand the coil behavior under the most likely SPARC TQ scenario. In addition to exploring how the MGI-induced TQ will affect the REMC effectiveness in the CQ (by way of current profile modification), we also consider how a closed-circuit REMC might affect the MGI-induced TQ, and we introduce further variation to the treatment of the region outside the SPARC limiter, showing that, as with the boundary conditions, this treatment can have significant consequences for the internal current profile evolution.

Like the simulations in Ref. 44, the simulations in this section are performed in a smoothed version of the SPARC vacuum vessel geometry [Fig. 9(a)] and have resistive wall boundary conditions for the $n>0$ components of the field, but an ideal wall for the $n=0$ component, inhibiting vertical instability. (As discussed previously, SPARC is expected to be passively stable to vertical motion with a 50 ms decay time,^49,50 but in the event of vertical motion in the CQ phase, the interaction with the REMC was found by Battey⁴⁹ to initially increase due to the motion of the plasma toward the horizontal coil segments, with a maximum coil coupling at 35 cm off the midplane). The magnetic perturbations associated with the increasing coil current are produced by the application of time-dependent normal magnetic field boundary conditions that produce the same perturbations at the last closed flux surface as those produced by the 3D coil. The resistive wall does not screen these external perturbations (the wall is transparent to the external perturbations but resistive with respect to the plasma response). The time-dependent amplitude is not pre-prescribed but grows in response to the decay of the plasma current. More details on these aspects of the model can be found in Ref. 44.

In addition to the MGI source, another new feature of these REMC simulations compared to those in Ref. 44 is the treatment of the region between the SPARC limiter and the vacuum vessel. This region (referred to as the ex-limiter region, as opposed to the plasma region) is treated as more vacuum-like in the simulation by the application of different physical parameters outside of the limiter than inside, such as larger resistivity, with a sharp tanh transition between them. As will be discussed in Sec. III B [and indicated in Figs. 9(b) and 9(c)], different parameters are enhanced in the ex-limiter region for different simulations. (This differs in several respects from the simulations by Sovinec in Ref. 52, which employs a different resistive wall model; in that case, separate equations are solved in explicitly separate plasma and vacuum regions, and the vacuum region is exterior to the resistive wall but interior to an outer-most ideal wall). In the divertor regions, the plasma-like parameters extend all the way to the simulation boundary. This enables open-field line currents, such as halo-currents, to flow to the boundary, reflecting the reality that such currents typically follow a return path through the vessel.

The REMC fields are applied at the boundary in the same way as in Ref. 44, with the REMC amplitude increase proportional to the decaying plasma current. In reality, the current in the REMC is driven by the loop voltage, and the change in plasma current is used as a proxy for this. Unlike the self-similar current decay prescribed in the ThinCurr code,⁴⁹ in which the coil current response is calculated, the current profile can redistribute in NIMROD, altering the internal inductance. Still, mapping coil current as a function of plasma current from ThinCurr is expected to be a closer approximation in NIMROD than simply using the ThinCurr coil current vs time, since the current decay time itself is not constrained in NIMROD to precisely match the ThinCurr calculation. In some cases, because it may be desirable (or necessary) to eliminate the complication of a passive switch in future devices, the REMC is permitted to carry current from the start of the simulation, mimicking a closed-circuit coil and allowing small perturbations from the coil even in the TQ. In other simulations, the coil perturbation is kept at zero until the time of the $Ip$-spike, behaving like a switch-activated coil for which a certain threshold voltage must be reached in order for the coil to carry current. All of the simulations in Sec. II D have the 6-valve non-resonant source distribution and the same 20% of deuterium in the source approximation as the high-resolution simulations in the previous section.

In Fig. 10, we compare the fluctuating field amplitudes for three MGI simulations with the REMC fields (a) off, (b) on only beginning at the $Ip$-spike, and (c) on throughout the simulation. In all simulations, we see that the dominant perturbation at early times is the n = 3 mode that is directly driven by the MGI source, and that a large n = 1 mode eventually becomes dominant. The MGI-only simulation is run for direct comparison, since the simulations in Sec. III A had a different boundary shape. Comparing this to the case with the REMC applied only after the TQ, we see that the large n = 1 mode (along with higher-n modes) initially begins to decay and can be seen to begin growing again due to the applied REMC fields beginning at around 4 ms. When the REMC is active from the start of the simulation, the n = 1 perturbations are evident in the pre-TQ phase, but smaller in amplitude than the n = 3 perturbations from the gas jet due to the relatively small changes in plasma current in this phase. However, these small n = 1 perturbations lead to an earlier growth and saturation of the large unstable n = 1 mode.

The effect of the closed-circuit REMC on the TQ is also seen in Fig. 11, in which the earlier onset of the TQ is again evident by the earlier peak in the radiated power, but with very similar amplitude. Differences also appear in the radiation toroidal peaking and the radiated energy fraction, with the latter being more pronounced. Both simulations have an increase in the radiation TPF at the time of the radiation peak, but a more pronounced brief peak in the TPF coincident with the radiation peak is evident for the closed-circuit case. In part, it is more pronounced because the TPF is lower before and after the radiation flash. More notably, the radiation fraction is higher for the closed-circuit case, rising from 95% to 99% by the end of the TQ when ionization energy is accounted for.

In Figs. 9(b) and 9(c), we see a comparison of the radiated power and plasma current for cases in which the parameters in the ex-limiter region are varied. In one case, the ex-limiter region is treated identically to the plasma region (“plasma-like”) although it is important to note that the resistivity in every case is non-uniform and varies at $T−3/2$⁠, where the ex-limiter region is typically clamped at a minimum temperature of 4 eV (⁠ $η=4.24×10−5$ Ohm-m for $Zeff=1$⁠). In the second case, the resistivity is enhanced by an extra factor of 1000 in the ex-limiter region. In the final case, the mass density (but not the number density) in the ex-limiter region is also multiplied by that factor (again there is no true vacuum region but a finite background plasma density evolving according to the same MHD equations throughout). This increased mass density serves to inhibit flows in the ex-limiter region, which could also be accomplished with increased viscosity. (Note that if charged particles flow across the boundary between these regions, momentum conservation can be broken in the enhanced mass density case; in practice, there is very little such flow). In each of these cases, the REMC becomes active at the time of the $Ip$-spike. The radiated power peak at the TQ exhibits different amplitudes and durations as the ex-limiter parameters change, with the largest amplitude being the plasma-like case. When both resistivity and mass density are enhanced, two peaks appear. This occurs because the enhanced mass density forestalls the growth of the large 1/1 mode, and a collapse of the temperature at mid-radius precipitated by higher-n modes triggers the first peak. This delay in the 1/1 onset is also responsible for the later $Ip$-spike in that case.

Even in the absence of large difference in the global current evolution, the internal current profile evolution is affected in a consequential way by the changes in the ex-limiter region parameters. A comparison of the q-profile as a function of time and normalized poloidal flux, $ψN$⁠, for the same three cases that appear in Fig. 9 as well as for a closed-circuit REMC case is seen in Fig. 12. The primary differences seen when the ex-limiter parameters are changed are, perhaps surprisingly, close to the core. At mid-radius, the safety factor rises sharply at the time of the $Ip$-spike, with the q = 2 surface shifting significantly inward. With the plasma-like ex-limiter region (a), the inward shift of q = 2 pauses at $ψN=0.25$ and only after another 1.5 ms does the on-axis q rise above 2. When the resistivity is enhanced in the ex-limiter region (b), this dynamic changes so that the $q0$ almost immediately rises above 2 at the time of the $Ip$-spike. However, when the mass density is also enhanced (c), the dynamic more closely resembles the original plasma-like ex-limiter case (a). Finally, the case with the closed-circuit REMC (d), which is otherwise the same as case (a), shows essentially the same behavior as (a) except for the slight change in timing due to the earlier TQ as seen previously. As, in general, the behavior of nonlinear simulations cannot always be easily interpreted, the reason for the effects of the edge parameters on the core current profile is not immediately clear, but the implication of these variations will be discussed further in Sec. IV.

The importance of maintaining $q0<2$ to maintain resonance of the plasma with the n = 1 REMC has been shown previously.^42,44 The fields are fully stochastic during the TQ phase in every case. Once $q0>2$ is reached, two islands typically appear, which subsequently coalesce into a core of good flux surfaces. A comparison of re-formed flux surfaces at 4.6 ms for two cases with the REMC applied at the time of the $Ip$-spike, but with different q-profile evolution, is seen in the upper and lower left of Fig. 13. In the lower case, flux surface reformation begins much earlier. In the upper case, flux surface reformation is just beginning now that $q0>2$⁠. We also compare this case with an identical simulation having no REMC, to see that, up to this time, the REMC has prevented the formation of a core of good flux surfaces that otherwise appears.

Two primary disruption mitigation systems designed for SPARC, the MGI system and the REMC, have been modeled with the NIMROD code. MGI into SPARC was modeled in four configurations with four or six valves, and all were found to produce high radiated power fractions and low radiation toroidal peaking. Based on a variety of input, including the NIMROD model, a six-valve configuration with three evenly spaced upper valves clocked $60°$ from three evenly spaced lower valves was selected for SPARC. High-resolution modeling of this configuration predicts a 2.5 ms pre-TQ with 5–15 GW of radiation, followed by a short TQ radiation flash of 65 and 20 GW of radiation in the early CQ. The radiation at the time of the TQ flash is localized to a single flux tube at mid-radius with a predominantly n = 1 toroidal variation and a TPF of less than 1.5 (and dominant wavelengths between 16 and 35 eV, see  the Appendix). At the end of the TQ, the radiation fraction is close to 99%, although this tends to dip during the CQ phase when the large 1/1 mode destroys confinement to such an extent that some thermal energy is conducted to the boundary.

In nonlinear modeling, the possible effects of varying simulation parameters such as diffusion coefficients can best be addressed with additional simulations, but generally speaking, the most significant approximations in the simulations herein are the high value of plasma viscosity and the constant value of parallel thermal conduction. Both of these approximations contribute significantly to convergence and numerical stability. Reduced viscosity could lead to both faster parallel transport and radial mixing of the impurities, which could result in radiation that is toroidally less peaked and more core localized, and therefore more favorable overall. Use of a large, constant value of parallel thermal conduction overestimates transport in colder regions of the plasma relative to the strong $T5/2$ dependence in the Braginskii model. Since parallel thermal conduction ultimately drives the parallel impurity spreading via the parallel pressure gradient, this would have an effect in the opposing direction to the enhanced viscosity, tending to underestimate toroidal radiation peaking. On the other hand, reduced parallel conduction in the cooler edge could limit the thermal energy conducted to the wall, particularly in the early CQ phase when the n = 1 mode peaks. Thus, the approximations are expected to have both favorable and unfavorable implications for mitigation metrics. Encouragingly, the SPARC MGI simulations performed with M3D-C1, using a different set of approximations, find reasonable agreement with NIMROD both with respect to the peak radiated power and the high radiation fraction.

The SPARC REMC has been previously modeled separately from the MGI system.^41,42,44 This modeling has predicted that the coil can fully suppress the RE current under some scenarios,⁴¹ but the resonance between the n = 1 coil and the plasma depends critically on maintaining $q0<2$⁠, which, in turn, can be affected by the TQ evolution and boundary conditions.^42,44 Here, we model the REMC under the most likely TQ scenario, which is triggered by the SPARC MGI system. The effects of a closed-circuit coil on the TQ are studied, and the n = 1 perturbation is found to trigger an earlier TQ onset and produce modest changes in the radiation fraction and TPF. We also vary the treatment of the ex-limiter region parameters, which, like the boundary conditions, are found to be important for the q-profile evolution. In some cases, on-axis q remains below 2 for 1.5 ms after the $Ip$-spike, and in other cases, it rises above 2 almost immediately. With a calculated CQ rate of 1.75 MA/ms, even an interval of 1.5 ms in which the coil effectively deconfines REs could result in a $103$ reduction in the overall avalanche gain. Still, a coil that remains resonant with the plasma throughout the CQ is more desirable. NIMROD is designed to accurately capture plasma magnetohydrodynamics, but better models for the interaction with what surrounds the plasma including limiters and conducting structures may be essential to making predictions with higher confidence.

This work was supported by Commonwealth Fusion Systems. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a Department of Energy Office of Science User Facility using NERSC Award FES-ERCAPm3195. Part of the data analysis was performed using the OMFIT integrated modeling framework.⁵⁴ 

The authors have no conflicts to disclose.

V. A. Izzo: Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). B. Stein-Lubrano: Investigation (supporting). A. Battey: Investigation (equal). R. Sweeney: Conceptualization (supporting); Investigation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). C. Hansen: Investigation (supporting); Writing – review & editing (supporting). R. A. Tinguely: Conceptualization (supporting); Investigation (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The particular wavelengths at which the peak radiation occurs will affect how much is reflected by plasma facing materials, therefore the likelihood of melting. Therefore, additional analysis was done for one simulation at the time of the TQ radiation flash to estimate the radiation spectrum. We estimate the individual line intensities at the time of maximum radiation for the 6-valve non-resonant (high-resolution) MGI simulation as follows. The version of NIMROD used for these simulations employs a coronal non-equilibrium radiation model, as in Refs. 27 and 31. (A more recent version, from which results have not been published, incorporates ADAS collisional-radiative rates.) In this model, charge state equilibrium is not assumed, and so the density of each Ne charge state is stored separately and is evolved by its own continuity equation. However, with this model, the intensity of line radiation can only be separated by charge state and not into individual line intensities. However, the data which NIMROD stores can be post-processed to obtain an estimate of individual line intensities using additional data files from OPEN-ADAS⁵³ that contain Photon Emissivity Coefficients (PEC) as a function of electron density and temperature for each Ne charge state. For instance, the OPEN-ADAS file $pec96#ne_pjr#ne0.dat$ is used for $Ne0$⁠, and the equivalent files for $Ne1+−Ne9+$ are used. Based on the individual charge state densities and electron density and temperatures at the time of maximum radiation, these data tables are evaluated at each point on the NIMROD 3D grid, and then, each individual line intensity is integrated over the volume to obtain a radiated power for that particular line. To save computational time, the volume integration is only performed for those lines with a maximum emissivity above a threshold that qualifies them as non-negligible. Of the 399 lines evaluated from the ten PEC files, 96 were found to have non-negligible radiated power, with a maximum value of 11.7 GW and a minimum of 0.40 MW. Radiation from these 96 lines sums to 70.1 GW. As a consistency check, this is reasonably close to the NIMROD radiated power of 65 GW at this time, given that a completely different set of radiation data are used, and the 399 data tables are evaluated with simple nearest-neighbor interpolation at each NIMROD grid point. In Table III, the radiated power for the largest 13 lines (those over 1 GW) is listed. A log –log scatterplot of all non-negligible lines is shown in Fig. 14.

The result of this analysis is that over 75% of the radiation occurs between 16 and 35 eV (35–77 nm). This result is close to the expected peak in Ne radiation based on the coronal equilibrium model,⁴⁸ which occurs near 30 eV.