The impact of collisionality, FLR, and parallel closure effects on instabilities in the tokamak pedestal: Numerical studies with the NIMROD code

Understanding the complex physics of instabilities in the tokamak pedestal region is essential for ensuring that tokamaks achieve their highest performance. Computation plays a critical role in furthering the knowledge of edge dynamics, including advances in the understanding of type-I edge-localized modes (ELMs).¹ Nonlinear computation of ELM dynamics remains challenging, but is becoming tractable with modern computing resources.^2–8 ELM modeling can be roughly broken into three regimes: the early linear regime, the nonlinear regime, and an energy loss regime. Two important aspects of the linear stage are the stability boundary within an operational space, e.g., pedestal current and pressure gradient, and the linear growth rate for each toroidal-mode number, i.e., the growth-rate or mode spectrum. The dominant instability drive that sets the linear-stability boundary is from the peeling-ballooning mode (PBM)^9–11 based on the ideal-MHD model. The linear spectrum is critical to the next phase, the nonlinear regime as it seeds the dynamics. Non-ideal effects can greatly impact the growth-rate spectrum and are thus important in determining the nonlinear ELM dynamics. In the final regime, high-fidelity simulation of the energy loss requires modeling of the fast conductive and convective transport, and accurate coupling to neutrals and the wall. State-of-the-art computations are beginning to push the boundary of current capabilities to modeling that captures accurate dynamics with these additional effects that are relevant during energy loss.^7,8,12

As the subsequent nonlinear evolution and energy loss in simulations are greatly impacted by the amplitudes of initial linear perturbations, a large source of uncertainty arises from the linear phase of ELM evolution. Given the difficulty (or impossibility) of measuring the small perturbations that seed an ELM, a reasonable model is to seed it from unstable PBMs that grow from small-wave perturbations. Although the stability boundary is well described by ideal MHD, details of the mode spectrum depend upon specific non-ideal effects. In particular, resistive effects may destabilize and first-order finite-Larmour-radius (FLR) effects may stabilize the modes. These latter effects, part of which is drift stabilization, require two-fluid modeling with advanced FLR closures. It is important to verify that codes accurately capture these effects within regimes of interest. In this paper, we study these effects on the linear modes with the extended-MHD NIMROD code.^13,14

Type-I ELM stability boundaries are dominated by drives described by ideal MHD as the ideal-MHD operator is stiff, i.e., near the ideal-MHD stability boundary, small changes in pressure and current-density gradients lead to large changes in the growth rate due to the ideal-MHD terms becoming large. Thus, only small profile changes are required to cross the boundary, and non-ideal modifications to the growth rate become secondary when describing the boundary location within operational phase space. This ideal-MHD stiffness is a general property that impacts the behavior of a wide class of plasma-instability phenomena. For example, in addition to ELM stability, it also describes the dynamics of disruptions caused by core-mode activity.^15,16 However, type-I ELMs are more complicated than core-mode disruptions given the broad toroidal-mode spectrum and associated nonlinear coupling. Similar to the studies of core-mode disruptions, comprehensive understanding of ELM onset requires evolving the equilibrium on the transport time scale through the instability boundary with a verified code that includes non-ideal terms to calculate an accurate initial linear-mode spectrum.

MHD codes that study ELM instabilities can be placed into two broad categories: (1) extended-MHD codes (e.g., NIMROD and M3D-C1¹⁷) and (2) reduced-MHD-based codes (e.g., BOUT++¹⁸ and JOREK¹⁹). The definition of extended MHD used in this work includes the full ideal-MHD force operator (see Ref. 20 and references therein). Using the full operator retains the fast compressional-Alfvén waves that enforce force balance and thus are critical to extended-MHD transport calculations with sources and sinks (e.g., Ref. 21). Relative to reduced modeling, the extended-MHD approach also has advantages in avoiding low-β approximations and a straight-forward implementation of the FLR closures as discussed in Section III. However, the inclusion of the full force operator comes at a cost of increased numerical difficulties when solving the equations associated with the fast waves.^22–26 Verification for all regimes is important, and previous work²⁷ verified the ideal-MHD terms for a non-diverted, peeling-ballooning unstable test case. As one expects that transport calculations across the ELM stability boundary coupled with an instability model that captures FLR drift stabilization significantly impacts dynamics, emphasis on verifying linear physics of extended-MHD codes thus naturally follows.

In this work, a verification of the full extended-MHD, an initial-value NIMROD code with the ideal-MHD ELITE code^10,11,28 is performed using the Meudas-1-case comparison.²⁹ This is a diverted tokamak case that is peeling-ballooning unstable. As we later discuss, the inclusion of the x-point significantly modifies the mode structure, making the comparison in full magnetic geometry an important additional verification. To assess the impact of the two-fluid, FLR effects, we compare our results with drift stabilization to a calculation of drift stabilization on the resistive PBM from Ref. 30. Although the approximations made during the analysis of first-order FLR PBMs in Ref. 30 preclude quantitative comparison, we make a qualitative comparison of our first-order FLR computations to analytics.

The paper proceeds as follows: In Sec. II, we summarize the benchmark and give an overview of our methodology. The exhaustive details of the methods used in this comparison, in particular, the mechanism whereby a vacuum region is effectively included through the extended-MHD modeling, are previously available in Refs. 27 and 31 and thus our description of this aspect is brief. In Sec. III, the model complexity is gradually increased from the ideal-like model through a resistive-MHD model with experimentally relevant profiles to a fully extended-MHD model which includes both parallel-closure and FLR-drift effects. The resulting impact on the growth-rate spectrum is analyzed. Comparisons of the computed drift-stabilized growth rates are made to analytic treatments of the PBM in Sec. IV. In Sec. V, the density and the temperature are varied at fixed plasma β (thus, the ion gyroradius, ρ_i, is modified) in order to ascertain how the FLR effects are impacted. In this study, the bootstrap current is fixed (the zeroth-order effect from variations of density modifies the bootstrap current) and the impact of FLR drift effects is isolated. Sec. VI summarizes this work and places it within the broader context of modeling of ELM relaxation events in the tokamak edge with transport effects.

We begin with a study of a high resolution, lower-single-null, JT-60U-like equilibrium (“Meudas-1”), which was originally employed in a benchmark of the MARG2D and ELITE codes,³² including a close approach to the X-point.²⁹ This extends the previous benchmarks²⁷ of ELITE and NIMROD as it includes diverted magnetic topology and a higher edge safety factor (⁠ $q 95 = 6.74$⁠, the safety factor at 95% of the normalized poloidal flux) that leads to increased resolution requirements. An ideal-MHD limit is achieved in NIMROD by using flat density and resistivity profiles inside the last closed flux surface (LCFS) with small resistivity (⁠ $S = 10 8$⁠), where S is the Lundquist number (⁠ $S = τ R / τ A$⁠), τ_A is the Alfvén time (⁠ $τ A = R o / v A$⁠), v_A is the Alfvén velocity (⁠ $B / m i n i μ 0$⁠), τ_R is the resistive diffusion time (⁠ $τ R = R o 2 μ 0 / η$⁠), $R o = 2.936 m$ is the radius of the magnetic axis, η is the electrical resistivity, μ₀ is the permeability of free space, $m α$ is a species mass (the α subscript denotes ions or electrons in this work), and $n α$ is a species density. The deuteron mass (m_i = 3.34 × 10⁻²⁷ kg) is used. In order to reproduce the vacuum response model outside the LCFS that is used by ELITE, a low density (0.01 of the core density) and a high resistivity (10⁷ times the core resistivity) are prescribed beyond the LCFS (more details on these approximations are in Ref. 27).

The normalized growth rates (⁠ $γ τ A$ where the linearized mode grows as $exp [ γ t ]$⁠) vs. toroidal mode number (⁠ $n ϕ$⁠) from NIMROD and ELITE are plotted in Fig. 1. There is good agreement between the codes except for $n ϕ = 4$ where there is a 27% relative difference. All other cases have a relative difference of less than 8% with typical differences of 5%. The NIMROD convergence in terms of the maximum polynomial order of the spectral elements is shown in Fig. 2. Convergence is most challenging at high wavenumbers where the resolution requirements are most stringent (the poloidal mesh is composed of 72 × 512 spectral elements). These cases converge from the unstable side where the growth rate decreases with enhanced resolution. Thus, the excellent agreement between NIMROD and ELITE at high $n ϕ$ in Fig. 1 may be spurious and indicate that slightly more resolution is required for $n ϕ > 25$⁠; however, the shown growth rates are likely within 5% of their converged values. Studying nearly ideal cases with the extended MHD codes such as NIMROD is challenging given the vanishingly small dissipation operators, and convergence is achieved more quickly with the additional non-ideal terms in the extended-MHD equations, as in the cases in Sec. III.

Figure 3 shows a poloidal cross section of the magnetic (B_R) eigenmode. The mode develops an “interference-pattern” structure near the X-point when inboard and outboard finger-like structures overlap. The finite-element-mesh nodes are superimposed atop the smallest-scale sub-figure. As established in Fig. 2, this simulation is spatially and temporally converged. The high resolution required to resolve these high-q₉₅, diverted cases motivated the development of memory-scaling improvements in the NIMROD code.

Initial-value extended-MHD computations, while more computationally expensive than ideal-MHD calculations with ELITE, are able to include additional effects that are more representative of experimental conditions:

• Instead of a magneto-static vacuum outside the LCFS, NIMROD includes a cold-plasma region.

• Finite-dissipation effects such as resistivity, viscosity, and thermal conduction.

• Density and temperature profiles that prescribe the associated drifts and resistivity profile.

• Parallel-closure effects that represent large diffusivities oriented along the magnetic field.

• First-order FLR-closure and two-fluid effects.

The authors in Ref. 31 investigated the first three of these effects with the Meudas-1 case. We further this work by including the last two effects. The first three effects are also re-investigated and different modifications to the growth rates are found. These differences are expected as slightly different profiles for density and temperature are used in order to avoid spurious electrostatic modes with the two-fluid, first-order FLR model. The authors in Ref. 31 avoided these modes by employing only the single-fluid resistive-MHD model.

In calculations beyond the ideal model, the density profile is $n e ( ψ n ) = n e 0 ( p ( ψ n ) / p 0 ) 0.3$ where $n e 0 = 6 × 10 19 m − 3$ and ψ_n is the normalized poloidal flux. Here, p is the pressure as prescribed by the ideal gas law, $p = n e T e + n i T i$⁠, where $T α$ is the species' temperature. This choice does not match Ref. 31, but effectively avoids spurious electrostatic modes with the two-fluid, first-order FLR model. These spurious modes are discussed in more detail in Sec. V. The pressure-profile gradient is specified by the Meudas-1 case, and setting the edge temperature to 200 eV leads to a core temperature of 5806 eV. The edge temperature is modified independent of the MHD-stability (determined by $p ′$⁠) by the transformation $p ( ψ n ) → p ( ψ n ) + p c$ where p_c is a constant. Without scrape-off-layer (SOL) profiles of density and temperature, which are not included in standard reconstructions, this choice of the edge temperature also specifies the temperature outside the LCFS. The SOL profiles are required to get both the correct profiles outside the LCFS, which determines the vacuum response, and simultaneously set the correct resistivity at the mode resonant location when Spitzer resistivity is used. The non-ideal viscous, conductive, and particle diffusivities are set to a small value, one tenth of the resistivity at the LCFS.

Figure 4 plots the growth rates computed by NIMROD on the “Meudas-1” case with four different models: ideal MHD (ideal); resistive MHD with reconstructed density, temperature; Spitzer-resistivity profiles (1f); and the single-fluid model with parallel Braginskii closures (1f-Par) where the parallel viscosity is

and the parallel heat-flux vector is

and a two-fluid model with parallel closures that includes ion gyroviscosity

where W is the rate-of-strain tensor

a fully extended MHD Ohm's law

and separate temperature evolutions with cross-heat fluxes

(2f-Par-FLR). Here, v_i is the bulk-ion flow, $ν ∥ i$ is the ion parallel diffusivity, $χ ∥ α$ is the species' parallel diffusivity, $q α$ (e) is the species' (electron) electric charge, J is the current density, $b ̂$ is the magnetic unit vector (⁠ $b ̂ = B / B$⁠), and I is the identity tensor. An effective ion charge (Z) of three is assumed. Each model builds on the previously listed and includes all prior terms. Electron viscosity is not included.

When a density profile is included, the normalized growth rate of the mode remains constant with the single-fluid model when the Alfvén time is computed at the mode resonant surface. Our calculations are normalized by the Alfvén time as computed with the values at the magnetic axis. These values remain constant when a density profile is included such that the normalized mode growth rate is effectively increased.

The effect of the resistivity profile is more nuanced. Finite resistivity allows for reconnection within the model by relaxing the frozen-flux constraint. This permits resistive-ballooning modes that grow faster than their ideal counterparts. Alternatively, the resistive dissipation can stabilize the modes by acting as a dissipative term and modifying the response outside the LCFS from that of a vacuum to that of a plasma. Figure 5 of Ref. 31 shows that the plasma response is stabilizing relative to a vacuum region. Overall, these effects stabilize the ballooning modes at low- $n ϕ$ and destabilize the modes at high- $n ϕ$ as seen in Fig. 4 when comparing the single fluid and ideal normalized growth rates. The effect of including small particle, viscous, and thermal diffusivities on the resistive-MHD mode is small (not shown).

Including the Braginskii parallel closures with large coefficients (⁠ $χ ∥ e = 10 9 × η 0 / μ 0$ and $χ ∥ i = ν ∥ i = 10 8 × η 0 / μ 0$⁠) has little effect on the growth rates. When thermal conduction is large, the temperature quickly equilibrates along field lines to produce an isothermal response. Thus, NIMROD computations that show little effect from the parallel closures are consistent with the ELITE results that find the shape of the mode growth-rate spectrum is largely not sensitive to the value of the ratio of specific heats, Γ, as seen in Fig. 1. When the response to the compressible motion associated with the sound wave is eliminated (Γ = 0 in the ELITE computations), the growth rate is slightly enhanced relative to the adiabatic limit. If the growth rates are modified beyond this small effect, it would be an indication that changes in the energy equation can impact the MHD response, and thus that closure effects are significant. Although the short-mean-free path Braginskii-like closure^34,35 is used in our NIMROD computations, we note that other closures, such as long-mean-free path³⁶ or the general approach of solving drift-kinetic equation³⁷ for a closure, can be applied. However, for this case where the parallel-closure effects are negligible, different parallel closures are unlikely to significantly modify the results.

With the full two-fluid, FLR model, there is a stabilizing effect on the intermediate and high- $n ϕ$ modes as expected from analytic treatments.³⁰ A common reduced-MHD approximation for first-order FLR closures is to assume that the ion-gyroviscous force and the advection by the diamagnetic drift exactly cancel (this is colloquially known as the gyroviscous cancellation, see, e.g., Ref. 38). However, as discussed by Ramos in Ref. 39, “these cancellations are only partial and not very useful in practice for general magnetic geometries…”. This cancellation is only valid with a large, uniform guide field, without curvature terms (slab approximation), and in the low-β limit (see, e.g., Ref. 40). The use of the full gyroviscous operator is critical as in addition to the effects that appear as diamagnetic drifts, it also contains terms proportional to first-order FLR magnetic-curvature and grad-B drifts.⁴¹ Our modeling contains not just the first-order FLR drifts from ion gyroviscosity but also the drifts from the cross-heat fluxes that enter the equations on the same order.

Just as the large balance of the ideal-MHD forces can cause numerical difficulties with the ideal-MHD force operator, the partial gyroviscous cancellation leads to similar numerical pitfalls and verification is important. The implementation of these terms in NIMROD is verified in cylindrical and slab magnetic geometries against analytic calculations for tearing^14,40,41 and ion-temperature-gradient (ITG)⁴² modes. To test the effect of these terms against analytic theory in full tokamak magnetic geometry, we perform a qualitative comparison with analytic theory in Sec. IV. Full verification is precluded by the approximations made in the analytic theories.

The first-order FLR model is valid for the toroidal modes presented, although this is not necessarily the case for general edge mode modeling. Assuming $k ≃ ( q / r + 1 / R ) n ϕ$ and evaluating local quantities at the peak of the mode eigenfunction (⁠ $ψ n = 0.969$⁠) on the outboard midplane, the normalized ion-gyroradius (⁠ $ρ i = Γ m i T i / Z e B$⁠) as a function of toroidal mode is $k ρ i ≃ 0.0055 n ϕ$⁠. If the first-order FLR assumption is violated, more-complex full-ion-orbit kinetic modeling is required to simulate large-fluctuation ELM dynamics. For this case, computations produce approximately the same results with and without the cross-heat fluxes. This is consistent with drift analytics⁴⁰ that shows that the cross-heat flux is only significant at very low values of plasma β (here the $β = 2 μ 0 p / B 2$ is 0.1%). Perhaps, coincidentally for this case, the effects of the resistive destabilization and drift stabilization largely counteract one another. The drift stabilized growth rates are not less than the ideal calculation until $n ϕ ≳ 24$⁠.

To gain insight into our two-fluid computations, we compare with analytic descriptions of drift stabilization of the ballooning mode. The dispersion relation for drift stabilizations of the ideal, incompressible ballooning mode is

where γ_I is the ideal growth rate in the absence of drift effects and ω is the complex frequency (⁠ $ω = i γ$⁠).⁴³ The diamagnetic drift velocity is typically defined as the perpendicular-to-the-magnetic-field flow contribution from the pressure gradient within the context of a two-fluid or generalized Ohm's law. For a tokamak, neoclassical effects damp the poloidal diamagnetic flow contribution. Thus in our case, the remaining toroidal diamagnetic flow determines the diamagnetic frequency, $ω * i = k · v * i = n ϕ / ( n i q i ) ∂ p i / ∂ ψ$⁠, where k is the mode wavenumber. Figure 5 shows the comparison to this dispersion relation with full and reduced diamagnetic drift values. The value of $ω * i$ is computed at the radial location where the mode structure peaks. The normalized frequency from the diamagnetic drift velocity, $ω * i τ A$⁠, varies linearly with toroidal mode number as $0.00483 n ϕ$⁠. This “local” approximation, which neglects the radial variation of the $ω * i$ profile, is known to over-estimate the effect of drift stabilization.^44,45 We find that a reduction of the value of $ω * i$ (approximately by a factor of 2–4) gives the best agreement between Eq. (7) and the two-fluid NIMROD calculation. A similar comparison of ELITE to the two-fluid BOUT++ model on a different case finds the drift stabilization is over-predicted⁴⁵ by a comparable factor.

One of the limitations of Eq. (7) is that it does not include the effect of resistive destabilization on the high- $n ϕ$ modes, as is included in our two-fluid computations. The analytic treatment of Ref. 30 provides a dispersion relation (Eq. (43) of the reference; referred to as HRP Eq. (43) in this discussion) that includes the effect of resistive destabilization in addition to drift effects. In order to compare with this theoretical treatment, parameters are computed on the outboard midplane at the radial location where the mode structure peaks, consistent with the previous $ω * i$ calculation. This yields local values of the safety factor, q = 7.34, normalized magnetic shear, s = 8.8, and normalized pressure gradient, $α = 28.8$⁠, where the definitions of Ref. 30 are used. The $Δ B ′$ stability parameter is determined by using the ideal growth rate (⁠ $Γ = 5 / 3$⁠) and solving Eq. 35 of Ref. 30. $Δ B ′$ varies with toroidal mode number but falls in the range of −17.4 to −12.1.

Comparison of HRP Eq. (43) with $ω * = 0$ and the NIMROD resistive-MHD results differ by a factor approximately 2.5 at large toroidal mode numbers where the analytic calculation results in greater growth rates and there is no stabilization at low $n ϕ$⁠. As the analytic calculation is a local calculation, it does not capture the stabilizing influence on the global mode of finite resistivity, in particular, it misses the effect of modifying the response outside the LCFS from that of a vacuum to that of a resistive plasma. Our interest is largely to compare and contrast the predicted drift stabilization, and thus, we reduce the resistivity to 35% of the value used in the NIMROD computations when calculating the analytic values of HRP Eq. (43). This produces growth rates with $ω * = 0$ that roughly match the NIMROD resistive-MHD computations at high $n ϕ$ as shown in Fig. 6.

More significantly, Fig. 6 compares the results of our linear NIMROD computations with the drift-stabilized growth rates computed from Eq. (43) of Ref. 30. Again, the analytics of HRP Eq. (43) over-estimate the effect of drift stabilization where the reduction of $ω *$ that leads HRP, Eq. (43), to qualitatively agree with the two-fluid NIMROD computations is between 2 and 4 (dependent on $n ϕ$⁠). There is little effect on the low- $n ϕ$ modes and the ion drift-wave resonance predicted in Ref. 30 is not observed.

In order to ascertain the role of two-fluid effects further, we vary the density and temperature while scaling parameters to maintain constant β and S (other dissipation parameters are scaled relative to resistivity). In experimental discharges, the zeroth-order effect on the edge of modifying the density is to modify the plasma collisionality and thus the bootstrap current. This causes a transition from ballooning-like modes at high density to peeling-like modes at low density. Our computations use a fixed equilibrium current and thus are not sensitive to this effect. Instead, we isolate the effect of the modification of ρ_i (and the associated ion skin depth, d_i) within the context of our two-fluid, first-order FLR model. This contrasts with Ref. 46 that examines these effects in concert and Ref. 47 that considers only the effect of modifications to the current profile.

Figure 7 shows the linear growth rate from NIMROD computations for resistive-MHD and three two-fluid cases with varying densities (the core density is varied from $3 × 10 19 m − 3$ to $1.2 × 10 20 m − 3$ where $6 × 10 19 m − 3$ is the value used in the prior computations discussed in the paper). The effect of drift-stabilization is stronger at low densities (ρ_i scales as $1 / n i$ at constant β) as expected. The cases with $n e = 3 × 10 19 m − 3$ are not plotted for $n ϕ > 18$ as the most unstable mode is no longer the peeling-ballooning mode but rather a dominantly electro-static mode related to the ITG as studied in the context NIMROD two-fluid advance in Ref. 42. This region is excluded as examination of ITG mode dynamics is outside the scope of this work.

Our comparisons to the ideal PBM growth rates computed by the ELITE code and analytic descriptions of the resistive, two-fluid PBM provides another verification test of the NIMROD code as a means to study the tokamak edge. Quantitative agreement with ELITE is achieved for these cases, extending prior work²⁷ that did not include diverted magnetic topology.

Two-fluid and FLR-drift effects can both produce stabilizing (drift stabilization through the sheared motion of the ion and electron fluid responses) and destabilizing (species decoupling and enhanced growth through mediation by the more mobile electron fluid) effects (see, e.g., Refs. 40 and 41). At intermediate and large mode numbers, the inclusion of finite-resistivity is destabilizing, while the FLR-drift effects are stabilizing. These effects are largely offsetting, resulting in a toroidal-mode growth-rate spectrum with the full extended-MHD model that resembles the ideal calculation. This result is not general but rather case dependent. Consistent with ELITE calculations that vary the ratio of specific heats from the adiabatic through the isothermal limit, we find the addition of anisotropic thermal conduction produces only a small modification of the growth rate.

Relative to the NIMROD calculations, analytics³⁰ over-predict both the effects of resistive destabilization and FLR-drift stabilization for the case studied. Qualitative agreement with the analytic descriptions of the peeling-ballooning mode is found where high- $n ϕ$ modes are susceptible to resistive destabilization and drift stabilization in both the analytics and NIMROD computations. One of the limits of the analytic derivation is a “local” approximation which leads to the over prediction of FLR-drift stabilization.⁴⁴ Given the highly complex nature of both the configuration of the tokamak edge and the two-fluid model equations, quantitative agreement is not expected as the local approximation neglects the global configuration and profile effects.

Regarding our investigation on the impact of the size of the ion gyroradius within the context of NIMROD's full first-order FLR, two-fluid model, we find significant changes on the growth rates of high- $n ϕ$ modes at varying densities but constant β. However, the most unstable mode number is largely unchanged. Thus, we conclude the zeroth-order effect of modifications to the bootstrap current largely dominates the determination of the most unstable mode during collisionality scans, while lower densities (large ρ_i) increase the magnitude of the drift stabilization on the high- $n ϕ$ modes consistent with the discussion of Ref. 48.

We thank Carl Sovinec, Chris Hegna, and Nate Ferraro for discussions involving this paper and N. Aiba for providing the Meudas-1 equilibrium. This material is based on the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award Nos. DE-FC02-06ER54875 and DE-FG02-08ER54972. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.