Toward integrated multi-scale pedestal simulations including edge-localized-mode dynamics, evolution of edge-localized-mode cycles, and continuous fluctuations

The tokamak edge region encompasses the boundary layer between hot core plasma and cold plasma facing components (PFCs). Therefore, this edge region includes complex geometry across the magnetic separatrix and rich physics, including plasma and atomic physics, as well as material science. The tokamak edge region is also vital to magnetic fusion program, which sets key engineering constraints for fusion reactor designs on power handling capability of divertor target and PFCs, and global energy confinement.

Depending on the heating power, plasma shaping, and magnetic configurations, various operating regimes have been observed in tokamak experiments, such as low-confinement mode (L-mode), high-confinement mode (H-mode),¹ and I-mode.² The low-confinement L-mode discharges have relatively poor particle and energy confinement. The high-confinement H-mode discharges have both good particle and energy confinement with the formation of an edge pressure pedestal. However, this mode of operation has been typically accompanied by periodic bursts of ejected plasma across the magnetic separatrix into the scrape-off-layer (SOL) labeled edge-localized-modes (ELMs).³ I-mode is a confinement regime that appears distinct from both L-mode and H-mode, with good energy confinement like H-mode, but with L-mode-like poor particle/impurity confinement. Therefore, I-mode has low pressure pedestal gradient with Edge-Localized Modes (ELMs) typically absent. Using external controls, such as shaping, ion $∇B$ drift direction relative to the X-point, plasma heating and fueling, these self-organized states of edge plasma exhibit different edge plasma density and temperature profiles and different underlying turbulent fluctuations to regulate their transport. Here are a few of such examples: broadband plasma fluctuations in L-mode and pre-H-mode,⁴ weakly coherent-modes (WCMs) in I-mode,⁵ quasi-coherent mode in Enhanced $Dα$ (EDA) H-mode in C-Mod,⁶ Edge Coherent Mode (ECM) in long-pulse H-mode plasmas in EAST,⁷ bursting chirping mode (BCM) in Enhanced H-mode pedestals with lithium injection in DIII-D,^8,9 edge harmonic oscillation (EHO) in quiescent H-mode (QH-mode) in DIII-D,¹⁰ and inter-ELM fluctuations labeled quasi-coherent fluctuations (QCFs) in ELMy H-mode.^11–13 Through further edge fluctuation spectrum control via E × B flow shear, a transition can be made from EHO to a regime with broadband turbulence, where the increased broad turbulence can reduce the pedestal pressure gradient, allowing the development of a broader and thus higher transport barrier in QH mode without ELMs.¹⁴ 

A detailed understanding and predictive capability for the turbulence and ELM instability serves not only as a projection of expected operating regimes but also as a guide for the development of regime access, ELM mitigation, and control capabilities. Over the last several decades, significant progress in the development of theoretical and comprehensive computational models of the turbulence and ELM has been made. BOUT (BOUndary Turbulence) was originally developed at LLNL in the late 1990s for modeling tokamak edge turbulence.^16–18 In collaboration with Univ. York, BOUT++ is a successor to BOUT and has been developed since late 2000s.¹⁹ Written in C++ to provide a friendly code interfacing environment, BOUT++ code is an object-oriented and a truly modular software framework for writing fluid/plasma simulations in curvilinear geometry on parallel computers. This open-source framework enables fast testing of numerical methods, quick development of a variety of advanced applications, and facilitating large team collaboration. A short list of related BOUT++ activities and presentations reported at this meeting includes the following topics. (1) A suite of two-fluid multiple-field models has been implemented in BOUT++ for all ELM regimes and fluid turbulence. (2) A suite of 3 + 1 Gyro-Landau-Fluid (GLF) models is also developed for pedestal kinetic turbulence and transport.

The high-fidelity BOUT++ two-fluid code suite has demonstrated significant recent progress toward integrated multi-scale simulations, including ELM dynamics, evolution of ELM cycles, and continuous fluctuations. Here is a list of principal results in this paper, along with comparisons with experimental observations. Section II gives an overview of an ELM crash dynamics in circular geometry. Nonlinear ELM simulations show three stages of an ELM event: (1) a linear growth phase; (2) a fast crash phase; and (3) a slow inward turbulence spreading phase. BOUT++ simulation results show a collisionality scaling of ELM energy losses consistent with the ITPA multi-tokamak database. The ITPA stands for International Tokamak Physics Activities. Nonlinear integrated multi-scale ELM simulations are described in Section III in circular geometry. A new coupling/splitting model has been developed to perform simulations of multi-scale ELM dynamics. BOUT++ simulation tracks five ELM cycles for 10 000 Alfvén times for small ELMs. To validate BOUT++ simulations against experiments and develop physics understanding of the fluctuation characteristics for tokamak operation regimes, both quasi-coherent fluctuations (QCFs) in ELMy H-modes and Weakly Coherent Modes (WCMs) in I-modes have been simulated in x-point geometry. The results are reported in Section IV. H-mode simulations predict that the H-mode profiles giving rise to QCFs are near marginal instability for ideal peeling-ballooning (P-B) modes, and the calculated particle diffusivity is smaller than the heat diffusivity. I-mode simulation results show that a strong instability exists at $n≥20$ for resistive ballooning modes and drift-Alfvén wave, and the calculated particle diffusivity is larger than the heat diffusivity, consistent with a key feature of the regime.

The Edge-localized Modes (ELMs) are quasi-periodic relaxations of the pedestal, resulting in a series of hot plasma eruptions, which potentially damage the ITER divertor plates and first walls. The image from a high-speed video in visible light of the MAST plasma obtained during ELMs shows ELM filamentary structures, which are parallel to the background magnetic field.¹⁵ Thomson scattering measurements from DIII-D show that the resulting ELM instability causes a rapid relaxation of the edge density and temperature profile. The accompanying divertor $Dα$ signal is a sequence of periodical pulses during ELMs, indicating the burst of heat flux flowing into the divertor.¹⁵ 

To study the physics of linear and nonlinear ELM dynamics, we create a set of JET-like circular equilibria using CORSICA²⁰ with the self-consistent variation of density and temperature profiles, while keeping the plasma cross-section shape, total stored energy, total plasma current, and total pressure profile fixed. The edge current is calculated from the Sauter bootstrap current model.²¹ We examine eight cases with increasing central density. As density decreases by a factor of 20, the neoclassical collisionality $ν*$ at peak gradient position decreases by a factor of 3000.

To simulate the ELM crash dynamics, we use the BOUT++ 3-field 2-fluid reduced MHD model, which evolves perturbed (1) pressure P, (2) vorticity ϖ, and (3) magnetic vector potential $A∥$⁠.²² The BOUT++ computation region in normalized poloidal flux ψ is $0.1≤ψ≤1.4$⁠, where ψ = 1 is the position of the magnetic separatrix. The complete description of the BOUT++ simulation model and results is given in our previous publication.²³ Here, we summarize the key results to illustrate the dynamics of the ELM crash and provide additional analysis for physics understanding of the ELM energy loss fraction vs collisionality scaling.

From the BOUT++ ELM simulations, the ELM energy loss fraction vs time is given in Fig. 4 of Ref. 23 and in Fig. 3 of Ref. 24, showing three stages of an ELM event: (1) a linear growth phase, (2) a fast crash phase, and (3) a slow inward turbulence spreading phase. In 3-field 2-fluid model, total energy loss shows a similar spreading as the 3 + 1 GLF model. In addition, in the 3 + 1 GLF model, the electron perturbation can provide additional spreading, which eventually dominates the total energy loss with a large conductive energy loss for cases with small ratio of the electron density scale length to the temperature scale length.

As the pedestal collisionality decreases, the bootstrap current increases, and the growth rate of the P-B mode decreases for high n but increases for low n (⁠$1≤n≤5$⁠). Therefore, the width of the growth rate spectrum Δ_n becomes narrower, and the peak growth shifts to lower n, as shown in Fig. 1. The 2D bi-spectrum in Fig. 2 shows that nonlinear mode coupling becomes stronger at high collisionality $ν*$⁠. Stronger nonlinear mode coupling at high $ν*$ leads to the lack of dominant filamentary structures and the reduced ELM energy loss.

ELM energy-loss fraction vs collisionality is shown in Fig. 4(b) of Ref. 23. BOUT++ simulations show collisionality scaling of ELM energy losses consistent with the ITPA multi-tokamak database, including DIII-D, JET, JT-60U with different shaping and heating methods.²⁵ We find that two fundamental factors determine if a single-mode amplitude can grow to a large magnitude to trigger an ELM: (1) linear growth rate spectrum and (2) nonlinear growth time. As the edge collisionality decreases, both linear and nonlinear physics set the ELM energy loss. Linearly, the dominant P-B mode shifts to lower n and the spectrum width of the linear growth rate decreases. Nonlinearly, the narrow mode spectrum leads to weak nonlinear phase scattering or weak nonlinear mode coupling as shown in Fig. 2, which leads to a long growth time and large ELMs.²⁶ The impact of the radial electric field E_r on peeling and ballooning modes is different. An increase in E_r significantly enhances the linear growth rate of low-n peeling modes, leading to large nonlinear ELM energy loss. An increase in E_r leads to large suppression of nonlinear ballooning amplitudes, but only weakly impacts their linear growth rates. The comprehensive analysis for the role of the radial electric field and its shear on the peeling-ballooning modes and its impact on ELM energy loss will be presented in a separate publication.²⁷ 

A new coupling/splitting model has been developed to perform simulations of multi-scale ELM dynamics. As a proof of principle, a set of three-field two-fluid model equations is split into axisymmetric components for background plasma profile evolution and fluctuation components for ELMs and the residual ELM-related turbulence. The axisymmetric component of the pressure is on the slow transport timescale, while the fluctuations are on a rapid timescale to changes in the profiles. The axisymmetric projection of nonlinear fluxes that are bilinear in fluctuating quantities, such as the energy flux, act as drive terms for the axisymmetric quantities that determine the profiles. No time evolving axisymmetric field components other than $〈p〉$ are included in the simulations reported in this paper. The equations are a set of coupled convection–diffusion equations, including (1) ELMs and the residual ELM-related turbulence, (2) flux-limited parallel transport, (3) micro-turbulence and neoclassical transport, and (4) sources and sinks. Some empirical models are introduced for items (3) and (4), which ensure that in the absence of the ELMs and the residual ELM-related turbulence, the axisymmetric pressure always returns to its initial profile. Because of the goal of this section is devoted to the development of multiscale algorithms, we do keep the choices of initial profiles and the length of time for recovery to test the schemes. The detailed description of the models is given in the  Appendix. Good agreement has been achieved with and without splitting for the same physics problem. With splitting, sources and sinks as well as items (2) and (3) have been added only to the axisymmetric component while keeping fluctuation components untouched.

As in Sec. II, we use JET-like circular equilibrium with the highest density case (⁠$n0=1×1020 m−3$⁠). Figure 3(a) shows the temporal evolution of ELM energy loss fraction, and Figure 3(b) shows the temporal evolution of the top pedestal pressure vs time. These two variables track each other for the pedestal pressure profile collapses and recover to a steep gradient during ELM cycles. The BOUT++ simulation tracks five ELM cycles for 10 000 Alfvén times, where the Alfvén time $τA=0.35 μs$⁠. The time traces from Divertor $Dα$ evolution in a DIII-D experimental discharge are shown in Figure 4(a), where the spikes in the Divertor $Dα$ signal indicate the burst of the heat flux flowing into the divertor after each ELM event. The time trace from electron pedestal pressure evolution is shown in Fig. 4(b). This is a typical DIII-D ELMy H-mode signature with large ELMs. Both BOUT++ simulations and DIII-D ELMy H-mode show the similar signature of a rapid relaxation of the pressure profile, even though the BOUT++ simulation is for small ELMs. Since we do not use the DIII-D experimental profiles and geometry, the comparisons are qualitative at best. Furthermore, the timescale is different between simulations and experiment because the simulations are for small ELMs and experiments for large ELMs.

The evolution of the pedestal pressure radial profiles is shown in Fig. 5. Figure 5(a) shows the pedestal pressure radial profile evolution during the recovery phase after an ELM crash from the BOUT++ simulation. The recovery time is 1225 τ_A. The simulations cannot recover to the initial dashed black curve due to the residual turbulent transport from ELMs; therefore, yellow curve is the steepest pressure profile for simulation recovery. Figure 5(b) shows the pedestal density radial profile evolution during recovery phase after an ELM crash from DIII-D fast Thomson scattering measurements. Both Figures 5(a) and 5(b) show that the pedestal pressure/density profile recovers to a steep gradient after an ELM event. However, the time scale does not match between simulations and experiments, because the simulations are for small ELMs and experiments for large ELMs. Therefore, comparisons of the pedestal and ELM characteristics with DIII-D experiments are again qualitative at best as we do not simulate the real DIII-D discharges or a shot with similar characteristics. Figure 5(c) shows the 2D radial-poloidal pedestal pressure profile during recovery phase at $t=500τA$⁠. There is a strong poloidal non-uniformity, even with the flux-limited Spitzer–Harm parallel heat diffusivity.

Depending on the heating power, plasma shaping, and magnetic configurations, edge plasmas can be in different operating regimes, such as low-confinement mode (L-mode), high-confinement mode (H-mode), and I-mode. In the different operating regimes, various coherent fluctuation structures have been observed in Alcator C-Mod, DIII-D, and other tokamaks. The edge turbulence dynamics is believed to play a key role for achieving different operating regimes. Here, we discuss those continuous fluctuations associated with ELMy H-mode and I-mode.

There is experimental evidence that quasi-coherent fluctuations (QCFs) are associated with the saturation of the pedestal between ELMs during ELMy H-mode discharges on C-Mod,¹¹ DIII-D,¹³ Asdex-Upgrade,²⁸ and JET.²⁹ QCFs have been observed as density and magnetic fluctuations. The QCFs are low-k, localized in the pedestal region with onset for a given temperature gradient, and track the evolution of the temperature gradient. QCFs are different from Quasi-Coherent-Modes (QCMs) observed in C-Mod ELM-free Enhanced $Dα$ H-modes in terms of frequency, wavenumber, and access conditions. On the other hand, the I-mode is characterized by changes in edge fluctuations compared to L-mode: (1) a Weakly Coherent Mode (WCM) and a broadband fluctuation suppression, (2) the WCM is localized close to the T_e pedestal (⁠$0.95≤r/a≤1$⁠), (3) the WCM and quasi-coherent mode (signature of the Enhanced $Dα$ H-modes) have similar radial localization and $kθ$⁠, but different frequency ramp up in different collisionality regimes.

In addition to the ELMs described in Sections II and III for a model JET-like circular geometry, fluctuations are simulated for the comparison of quasi-coherent fluctuations (QCFs) between ELMs in H-mode and weakly coherent modes (WCMs) during ELM-suppressed I-modes in full toroidal x-point geometry for BOUT++ validation and physics understanding. We found that (1) H-mode simulations predict that the QCFs are near marginal instability for ideal peeling-ballooning modes, and the calculated particle diffusivity is smaller than the calculated heat diffusivity; (2) I-mode simulation results show that a strong instability exists at $n≥20$⁠, for resistive ballooning modes and drift-Alfvén waves; the calculated particle diffusivity is larger than the calculated heat diffusivity.

For BOUT++ simulations, we use three-dimensional 6-field 2-fluid electromagnetic model, which evolves perturbed (1) ion density n_i, (2) temperature T_i, (3) parallel velocity $v∥,i$⁠, (4) electron temperature T_e, (5) vorticity ϖ, and (6) magnetic vector potential $A∥$⁠.³⁰ The simulations are based on a set of C-Mod and DIII-D experiments, including (1) fitted plasma profiles and equilibrium reconstruction utilizing experimental pressure constraints (kinetic equilibrium fitting code (EFIT)) and (2) BOUT++ computation region that spans the magnetic separatrix, including three distinct regions: the outer part of the closed flux region (Edge), the SOL, and the private flux region.

To simulate pedestal plasma turbulence and validate with the corresponding experimental measurements, the BOUT++ code uses realistic X-point magnetic and plasma profiles. The background magnetic field structure is obtained from an MHD equilibrium code (usually kinetic EFIT) for the shots simulated. The plasma profiles are obtained from measurements typically averaged over time. From the given magnetic geometry and plasma profiles corresponding to a specific experimental device and discharge, the simulation is initialized with a set of small random fluctuations. The fastest growing modes dominate the initial phase of the calculation, in which the fluctuations grow at an approximately exponential rate. After this initial linear phase, the density, temperature, electrostatic, and magnetic potential fluctuations evolve to a saturated state with many modes. From the saturated steady state, turbulence statistical properties can be extracted from the BOUT++ simulations using Fourier analysis to steady turbulence and validated with the various fluctuation measurements.

For simplicity, the profiles for the axisymmetric components have been held fixed in this type of simulation for comparison with experimental measurements. This assumption is reasonable, because the experimentally measured profiles are statistically similar in the saturation of continuous fluctuations in ELMy H-mode and I-mode phases. The equilibrium electric field is either from experimental measurement when available (such as in I-mode simulations) or determined from force balance with no net flow $Er0=(1/n0Zie)∇rPi0$ with ion pressure $Pi0$ (such as in H-mode simulations). We start the simulation using no net flow in $Er0$ term when flow measurements are not available. The turbulent zonal flow (ZF) has been set to be zero as well. The perturbed electric field is $Ẽr=(1/n0Zie)∇rp̃i$⁠. The impact of these assumptions is a subject of further investigation using the BOUT++ code. As is well-known, nonlinear interactions between turbulence and zonal flows (ZF) or geodesic-acoustic modes (GAM) have been extensively studied because of its important role during transitions of tokamak confinement regimes. Furthermore, it is also observed that both the WCM and the GAM are necessary for the L-to-I transition and co-exist in I-mode regime, and I-mode regime access is found to be sensitive to the GAM drive and damping.^31,32 The zonal magnetic field is also set to zero as it is negligibly small compared with the equilibrium magnetic field B₀. However, to simulate transition physics, in addition to matching the turbulence characteristics, extra efforts will be needed to ensure consistent sources and sinks added to match experimental radial profile evolution, zonal flows, and GAMs, which is beyond the scope of present paper. In the following, we first summarize C-Mod QCFs simulations and comparison with experimental measurements, then DIII-D QCFs, finally followed by the C-Mod I-mode.

Figure 6 shows the C-Mod plasma profiles in pedestal region for ELMy H-mode discharge 1120815027 at time 1075 ms (red) and I-mode discharge 1120907032 at time 1010 ms (blue). Here, (a) electron density, (b) electron temperature, (c) total pressure, (d) parallel current in pedestal region, calculated using the Sauter formula, and (e) Last-Closed-Flux-Surface (LCFS) contours. The ELMy H-mode discharge was performed with $BT=5.4 T$ and current $Ip=0.9 MA$⁠. In comparison with ELMy H-mode pedestals, I-mode has a lower pedestal pressure and current than H-mode due to lack of a particle barrier and, therefore, exhibits ELM-free operating properties. The I-mode regime is obtained typically with the ion $∇B$ drift away from the active X-point as shown in Fig. 6(e), in the “unfavorable” configuration for obtaining H-mode with significant higher power threshold to access H-mode.

BOUT++ simulations for C-Mod pedestal show that ELMy H-mode and I-mode exhibit different underlying instabilities. Linear simulations show that the I-mode profiles are unstable for resistive ballooning mode and drift-Alfvén wave as shown in Fig. 7(a), while the ELMy H-mode profiles are marginally unstable near the peeling-ballooning threshold as shown in Fig. 7(b). We also test the case without curvature shown in blue diamond and the growth rate of unstable modes drops by almost half. For I-mode simulation, the residual instability without curvature is the collisional electromagnetic drift-Alfvén waves, while for the ELMy H-mode, the residual instability without curvature is the peeling mode. The P-B stability of the H-mode and I-mode pedestal calculated in this way is consistent with a prior comparative analysis using the ELITE code.³³ 

Nonlinear BOUT++ simulations have been successfully performed for the evolution of quasi-coherent fluctuations. The spectrogram for quasi-coherent fluctuations of C-Mod ELMy H-mode is shown in Figure 8: (a) frequency vs radius from BOUT++ simulations; (b) edge $Hα$ trace indicates the ELMs; (c) frequency versus time from inter-ELM magnetic fluctuations measured close to LCFS using the Mirnov coils. Similar density fluctuation characteristics are also observed from reflectometer measurement. The evolution of the poloidal wavenumber spectrum for quasi-coherent-fluctuations of C-Mod ELMy H-mode is shown in Figure 9: (a) BOUT++ simulations; (b) two-point correlation measured using a double-head magnetic probe. BOUT++ simulations suggest that QCFs are localized in the pedestal peak pressure gradient region having a predominant frequency at $f≃300−400 kHz$ and poloidal wavenumber at $kθ≃0.7 cm−1$⁠, and propagate in the electron direction in the laboratory frame. The overall signatures of simulation results for QCFs show good agreement with C-Mod measurements. It is worth noting that the timescale of the evolution of the poloidal wavenumber spectrum is different between simulations and experiments. Because the plasma profiles are fixed in the simulations, once fluctuations evolve to a saturated state from a linear phase, the assumption is that the fluctuations will remain in the saturated state and their turbulent spectra also remain the same. Therefore, the mismatch of timescale between simulations and experiments is not relevant. Furthermore, because the plasma profiles are time-averaged between ELMs from experimental measurements, the comparison should be referred to the experimental spectrum during the inter-ELM cycles.

Similarly, we have performed linear and nonlinear BOUT++ simulations for quasi-coherent fluctuations (QCFs) on DIII-D Type-I ELMy H-mode discharges. The experiments were carried out on the DIII-D tokamak at fixed $BT=1.9 T$ for three different values of plasma currents in a lower single-null configuration with the ion $B×∇B$ drift direction toward the X-point.¹³ These experiments were targeted at capturing the pedestal recovery after an ELM crash using the recently upgraded Thomson scattering system with high temporal resolution using multiple lasers. Here, we conduct BOUT++ simulations for the high plasma current $Ip=1.6 MA$ case with the pedestal plasma profiles shown in Figure 10 for electron density, temperature, the total pressure, and parallel current. Linear analysis indicates (1) the pedestal is marginally unstable for ballooning modes and the collisional electromagnetic drift-Alfvén wave appears to be dominant. Nonlinear analysis indicates that the frequency of the dominant mode is about 80 kHz which is near to that of the observed QCFs as shown in Figure 11. The experimental measurements show both lower and higher frequencies of QCFs. The high frequency should be the 2nd harmonics of the low frequency, which is not observed in simulations.

Correlation of QCF amplitudes' evolution with the gradient of the pedestal electron temperature has also been observed on DIII-D,¹³ Alcator C-Mod,¹² and Asdex-Upgrade²⁸ for ELMy H-mode between ELMs. Shown in Fig. 12(a) is the evolution of the magnetic fluctuations and pedestal electron temperature gradients during ELM recovery on DIII-D. It is shown that QCF has an onset for a given critical temperature gradient $∇Te$⁠, and electron temperature gradient $∇Te$ and QCF amplitude track each other. When the gradient of the pedestal electron temperature increases to a critical value, it triggers the QCF to grow (marked as the onset of QCF in Fig. 12(a)). By the time when QCF grows to large amplitude at about 20 ms, we suspect that the large QCF transport balances the electron heating power from core plasma, and the gradient of the pedestal electron temperature stops growth and reaches the saturation. Based on these experimental observations, the QCF appears to limit the gradient of the pedestal electron temperature. To clearly illustrate the correspondence, we have performed linear BOUT++ simulations for a scan of electron temperature gradient $∇Te$ with the most unstable toroidal mode number n = 30, ion temperature, and total pressure profiles held fixed in C-Mod ELMy H-mode. However, the density profile is changed to offset the changes of electron temperature. The results are given in Fig. 12(b). It is clearly shown that increasing $∇Te$ at the pedestal leads to a strong increase in linear growth rate, also implying an increase of the QCFs intensity, while the linear growth is not sensitive to increasing $∇Ti$⁠, possibly due to the balance of two factors: increased ballooning drive offset by the ion diamagnetic stabilization. The insensitivity of the QCFs to the ion temperature gradient $∇Ti$ indicates that the ion temperature gradient $∇Ti$ does not track with QCF amplitude, since the ion temperature gradient $∇Ti$ is much more relaxed on DIII-D than the electron temperature gradient $∇Te$ as shown in Fig. 10. The nonlinear simulations for the scan show that when the gradient of the pedestal electron temperature $∇Te$ reduces from experimental value $∇Te, exp$⁠, the QCFs become less coherence. For example, the same coherence properties of QCFs remain for $∇Te=90%∇Te, exp$⁠, but the QCFs lose their coherence for $∇Te<70%∇Te, exp$⁠. Sensitivity studies of the nonlinear simulations for QCF amplitudes with respect to pedestal density, pedestal pressure, and parallel current profiles will be presented in a future publication.

The complete BOUT++ I-mode simulation results and extensive comparison with C-Mod measurements are given in a separate publication.³⁴ Here, we first summarize the main simulation results and show the comparison with experimental measurements. Then, to illustrate the key feature of I-mode pedestal with no particle barrier, we compare the calculated transport from nonlinear simulations for both weakly coherent modes (WCMs) in I-mode and quasi-coherent fluctuations (QCFs) in ELMy H-mode.

Figure 13 shows both mode number and frequency spectra. Figure 13(a) shows the contour plot of density fluctuation vs toroidal mode number n and time at the peak gradient position $ψ=0.98$ and outside midplane. It is clearly shown that the most unstable mode number n shifts from linear phase at n_max = 30 to nonlinear phase n = 20. The frequency spectrum of nonlinear BOUT++ simulation features a peak around $∼300 kHz$ for the mode number n = 20 that is similar to the reflectometer measurement at nearby position, as shown in Fig. 13(b).

From BOUT++ 6-field 2-fluid nonlinear simulations, we calculate the particle and heat diffusivities for electrons and ions separately (using the fluxes defined in Eqs. (11) and (12) of Ref. 35 by dividing the corresponding negative radial gradient); the results are given in Fig. 14 for I-mode and ELMy H-mode. For C-Mod I-mode, Fig. 14(a) shows the particle diffusivity (solid curve), electron thermal diffusivity (long-dashed), and ion thermal diffusivity (short-dashed). It is clearly shown that (1) the calculated χ_e agrees with the χ_eff from the experiment, and (2) the calculated particle diffusivity is larger than the heat diffusivity. The diffusivities for DIIID ELMy H-mode are shown in Fig. 14(b) with the same legends. In addition, classical and neoclassical particle diffusivities are shown in dotted-dash and dotted curves. For ELMy H-mode, the calculated particle diffusivity is smaller than thermal diffusivity (⁠$D≤χe$⁠) and is close to the neoclassical level. However, the calculated particle and heat diffusivities for C-Mod ELMy H-mode have similar amplitude.

The high-fidelity BOUT++ two-fluid code suite has demonstrated significant recent progress toward integrated multi-scale simulations that reproduce features of experiments, including ELM dynamics, evolution of ELM cycles, and continuous fluctuations.

Nonlinear ELM simulations in JET-like circular geometry show three stages of an ELM event: (1) a linear growth phase; (2) a fast crash phase; and (3) a slow inward turbulence spreading phase. BOUT++ simulation results show a collisionality scaling of ELM energy losses consistent with the ITPA multi-tokamak database. Two fundamental factors are found to determine ELM energy losses: (1) linear growth rate spectrum and (2) nonlinear growth time. As the edge collisionality decreases, both linear and nonlinear physics set the ELM energy loss. Linearly, the dominant P-B mode shifts to lower n, and the spectrum width of the linear growth rate decreases. Nonlinearly, a narrow mode spectrum leads to weak nonlinear mode coupling, which leads to a long growth time and large ELMs.

A new coupling/splitting model has been developed to perform simulations of multi-scale ELM dynamics in JET-like circular geometry. As an example, a set of three-field two-fluid model equations is split into axisymmetric components for background plasma profile evolution and fluctuation components for ELMs and the residual ELM-related turbulence. The axisymmetric component is a set of coupled convection–diffusion equations, including (1) ELMs and the residual ELM-related turbulence, (2) flux-limited parallel transport, (3) micro-turbulence and neoclassical transport, and (4) sources and sinks. Some empirical models are introduced for items (3) and (4), which ensure that in the absence of the ELMs and the residual ELM-related turbulence, the axisymmetric pressure always returns to its initial profile. With splitting, sources and sinks as well as items (2) and (3) have been added only to the axisymmetric component while keeping fluctuation components untouched. The simulation tracks five ELM cycles for 10 000 Alfvén times for a small ELM. The temporal evolution of the pedestal pressure is similar to that of experimental measurements for the pedestal pressure profile collapses and recovers to a steep gradient during ELM cycles.

To validate BOUT++ simulations and develop physics understanding of the fluctuation characteristics for different tokamak operation regimes, both QCFs in ELMy H-modes and WCMs in I-modes have been simulated for comparison with experiments using three-dimensional 6-field 2-fluid electromagnetic model. H-mode simulations predict that the QCFs are near marginal instability for ideal peeling-ballooning modes. BOUT++ simulations suggest that (1) QCFs are localized in the pedestal region having a predominant frequency at $f≃300−400 kHz$ and poloidal wavenumber at $kθ≃0.7 cm−1$⁠, and propagate in the electron diamagnetic direction in the laboratory frame. The overall signatures of simulation results for QCFs show good agreement with C-Mod and DIII-D measurements. (2) The ELMy H-mode pedestal is near the marginal instability threshold for ideal peeling-ballooning (P-B) modes for both C-Mod and DIII-D, while the collisional electromagnetic drift-Alfvén wave appears to be dominant for DIII-D. (3) Particle diffusivity is either smaller than the heat diffusivity for DIII-D or similar to the heat diffusivity for C-Mod. I-mode simulation results show that a strong instability exists at $n≥20$⁠, for resistive ballooning modes and drift-Alfvén waves. For nonlinear I-mode simulations, the frequency spectrum of nonlinear BOUT++ simulation features a peak around $∼300 kHz$ for the mode number n = 20 that is similar to the reflectometer measurement at nearby position. The calculated particle diffusivity is larger than the calculated heat diffusivity, $D≫χe$⁠, which is consistent with the key feature of I-mode pedestal with no particle barrier.

From this study, we demonstrate that the underlying linear instability drives and turbulent fluctuations have distinct characteristics for different tokamak operation regimes, which can be accessed by many external controls, such as shaping, ion $∇B$ drift direction relative to the X-point, plasma heating, and fueling. By looking forward to further development of physics understanding of these self-organized states of edge plasma, which exhibit different edge plasma density and temperature profiles and different underlying turbulent fluctuations to regulate their transport, we have to study the transition physics from one regime to another, including from QCFs to ELM crashes. To do so, in addition to matching the turbulence characteristics, extra efforts should be made in BOUT++ code development and validation to ensure consistent sources and sinks added to match the experimental radial profile evolution including toroidal flow, zonal flows, and GAMs as described in Sec. III.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. This material is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences. The authors wish to thank the support of the China Natural Science Foundation under Contract Nos. 11405215 and 11405217, 11505222, 11405214, and Y45ETY2306. National Magnetic Confinement Fusion Science Program of China under Contract Nos. 2014GB106001 and 2014GB106003, and 2015GB101000. The authors wish to thank B. Cohen, J. G. Chen, A. Dimits, B. Dudson, P. H. Diamond, M. E. Fenstermacher, T. Golfinopoulos, C. Holland, I. Joseph, M. Kim, S. S. Kim, D. F. Kong, G. Q. Li, C. H. Ma, J. F. Ma, B. H. Meyer, P. B. Snyder, C. K. Sun, T. F. Tang, M. V. Umansky, H. Wilson, Z. H. Wang, and P. W. Xi for useful discussions. LLNL-JRNL-680060.

As a proof of principle, we take a three-field model and split the equations into axisymmetric components for background plasma profile evolution and fluctuation components for ELMs and the residual ELM-related turbulence. A similar strategy can be generalized to other fluid and Gyro-Landau-Fluid equations.

The set of three-field model equations can be written as follows:

Here $∇‖F=B∂‖(F/B)$⁠, for any F, $∂‖=∂‖0+b̃·∇,b̃=B̃/B=∇‖A‖×b0/B,∂‖0=b0·∇,κ0=b0·∇b0$⁠.

Defining $p=〈p〉+p̃, ϖ=〈ϖ〉+ϖ̃$⁠, $A‖=〈A‖〉+A‖̃$⁠, and $ϕ=〈ϕ〉+ϕ̃$⁠, where $〈〉$ means the averaging over z (binormal coordinate) for toroidal mode number n = 0, Eqs. (A1)–(A3) can be separated into equations in the fluctuation and averaged (or axisymmetric) parts. The set of fluctuating equations is as follows:

Here, the total time derivative only contains the convective part of E × B velocity, $d/dt=∂/∂t+vE·∇$⁠. We also assume no equilibrium vorticity $ϖ0=0$⁠.

The set of axisymmetric (or averaged) equations is as follows:

Notice that the z-averaged (or axisymmetric) variables $〈p〉$ and $〈ϖ〉$ make no direct contribution to the averaged energy flux $Γ=〈p̃(∂ϕ̃/∂z)〉$⁠, where $〈p〉$ contributes to turbulence fluctuations through the flattening of the background density profile and $〈ϖ〉$ contributes through the velocity shear effect. Of the axisymmetric field components, only the time evolution of $〈p〉$ has been retained in the simulations reported in this paper.

The set of axisymmetric equations can include sources, sinks, and other transport processes for pressure, flow, and current. For example, with the sources, sinks, and other transport processes for pressure, the set of axisymmetric (or averaged) equations can be written as follows:

where $χψ$ and $χ∥$ are radial and parallel transport coefficients and only function of axisymmetric variables, such as $〈p〉$⁠. The parallel diffusion coefficient $χ∥$ is the flux-limited Spitzer–Harm parallel heat diffusivity. The radial diffusion coefficient $χψ$ can be neoclassical and anomalous. For a given equilibrium pressure profile $P0(ψ)$⁠, the radial diffusion coefficient $χψ(ψ)$ can be obtained by inverting the equilibrium flux-surface-averaged pressure equation with given sources, sinks, and boundary conditions. In a case with no ELM events, the ELM E × B convection $〈vE·∇p〉$ on the left-hand side in Eq. (A23) is zero, and the pressure equation has no coupling with other fluctuation equations. Therefore, in the steady state with no ELMs, Eq. (A23) maintains the initial pressure profile. Since the initial pressure profile is above the peeling-ballooning instability threshold, the ELM bursts, which leads to the initial pressure profile collapses and the residual ELM fluctuation leads to further turbulent spreading and quasi-linear flattening. To simulate the multi-ELM cycles, additional sources and possibly the suppression of the residual ELM-driven fluctuations are needed for pedestal rebuilds. Because the goal of this section is devoted to the development of multiscale algorithms, we do keep the choices of initial profiles and the length of time for recovery to test the schemes.

It is worth noting that the ELM fluxes $ΓP=〈vE·∇p〉$ contain both diagonal and off-diagonal turbulent transport driven by the radial gradients of axisymmetric vorticity $∇〈ϖ〉$ and axisymmetric pressure $∇〈p〉$⁠. It is the more general form than previous coupling strategies.^36,37

For a given ELM crash, there is a large perturbation of zonal pressure profile $〈p〉≠0$⁠, which flattens the total pressure profile (⁠$P=P0+〈p〉$⁠) and leads to an ELM crash. For a given neoclassical and anomalous diffusion coefficient $χψ$⁠, to simulate ELM cycles, we need to introduce an additional source p_source to compensate the additional ELM loss. For a simple case, we assume that the source takes the following simple form:

The axisymmetric pressure equation becomes

Here, J and g₁₁ are the coordinate Jacobian. If the zonal pressure perturbation is damped to zero via heating and residual core transport, the total pressure profile is restored to its initial pedestal profile; the ELM event repeats. Because we use no net flow assumption and there is force balance for both the equilibrium and the perturbation, the E × B flow shear will also be increased after the total pressure profile is restored to its initial pedestal profile.

In the nonlinear integrated multi-scale ELM simulations presented in Sec. III, the source coefficient $γsource=0.01$⁠.