Comparison of Saturation Rules Used for Gyrokinetic Quasilinear Transport Modeling

Abstract

Theory-based transport modeling has been widely successful and is built on the foundations of quasilinear theory. Specifically, the quasilinear expression of the flux can be used in combination with a saturation rule for the toroidal mode amplitude. Most transport models follow this approach. Saturation rules are heuristic and difficult to rigorously derive. We compare three common saturation rules using a fairly accurate quasilinear expression for the fluxes computed using local linear gyrokinetic simulation. We take plasma parameters from experimental H-mode profiles and magnetic equilibrium and include electrons, deuterium, and carbon species. We find that the various saturation rules provide qualitatively similar behavior. This may help to explain why the different theory-based transport models can all predict core tokamak profiles reasonably well. Comparisons with nonlinear local and global gyrokinetic simulations are discussed.

1. Introduction

Prediction of turbulent particle and energy transport is critical for improving the performance of a fusion reactor. Much progress has been made with reduced models in the core region, from the pedestal top inwards. Theory-based models, including the trapped gyro-Landau fluid (TGLF) model [1,2,3,4], the multi-mode model (MMM) [5,6,7,8], and the QuaLiKiz gyrokinetic transport model [9,10,11,12] are successful in predicting core density and temperature profiles over a range of tokamak plasma operating conditions. Quasilinear theory is used widely to compare with both experiment and nonlinear gyrokinetic simulation [13,14,15,16]. Typically, the quasilinear expression is taken for the flux and a heuristic saturation rule is invoked to obtain the mode amplitude, thereby determining the nonlinear flux. While the level of the fluxes obtained using this type of approach may not be accurate, the parametric dependence on wavelength and plasma parameters is often insightful. Even with the successes of the various theory-based transport models for predicting core density and temperature profiles, there remains a need for better understanding. For example, particle transport and associated density buildup is not well understood [17,18]. Additionally, High-Z impurities, e.g., Tungsten in ITER, do not fully ionize, and can produce significant radiative power loss if core concentrations are not well controlled [19,20].

Here, we further examine the quasilinear transport modeling approach. We compare it to local and global nonlinear gyrokinetic simulation in order to determine which can best model the governing equations with relatively few approximations. We directly compare three widely used saturation rules, two of which are derived from simple scaling arguments [9,14] and a third which has been shown to provide reasonable parameter dependence for fluxes [21]. While the comparisons we present are rudimentary, we are unaware of any such previous study of sensitivity to the saturation rule. In addition, we compare it with the TGLF model, which is specifically designed to agree with flux-tube gyrokinetic simulation from the CGYRO code [22,23]. The goal of this paper is to examine the sensitivity of the choice of saturation rule, which is a part of the theory that is least well understood. Generally, tokamaks operate within regimes that have relatively good confinement; hence, it is not unreasonable to assume that the turbulence is weak and made up of a number of active linear eigenmodes that interact due to weak nonlinear coupling. Linear calculations with gyrokinetic codes are routine and computationally fast. Linear fluxes can easily be obtained from gyrokinetic simulation. In Section 4, we use GENE linear calculations to obtain the quasilinear expression for the flux. Nonlinear simulations are much more compute-intensive. However, no information on the saturation level is available from linear calculations. Therefore, it is common to invoke a “saturation rule” that provides the saturation level of the turbulence as a function of the linear growth rate, wave number, and other parameters. The ability to derive a rigorous saturation rule is elusive. One reasonable approach is to obtain an empirical saturation rule using the scaling of nonlinear gyrokinetic simulation [2]. While the saturation rule is probably the weakest link with regard to rigor, the assumption of a quasilinear expression for the flux and how the quasilinear flux is calculated may be approximate.

In this paper, we discuss results from gyrokinetic simulation using the GENE and GEM codes. The GENE code is used for linear and nonlinear local calculations, including the calculation of the quasilinear expression for the flux [24,25]. We use GENE for linear calculations due to its high accuracy, good convergence properties, and comprehensive physics capability. GEM is an efficient tool for nonlinear global simulation due to its robust behavior over a wide range of parameters and its fast performance on parallel computing platforms [26]. In this study, we choose realistic plasma profiles and magnetic equilibrium from a conventional ELMy H-mode DIII-D case (162,940) just prior to the onset of an ELM [27]. We include electrons and two ion species, namely, deuterium (main) and carbon (impurity). The details are discussed in Section 2. We begin by discussing the plasma parameters for our study in Section 2. In Section 3, we investigate the linear properties of the selected profile. In Section 4, the quasilinear theory is described and the resulting turbulent transport is compared for three different saturation rules. In Section 5, we compare nonlinear fluxes from the GENE and GEM codes.

2. Tokamak Plasma Parameters and Assumptions

To compare the three saturation rules in quasilinear theory, we use DIII-D discharge 162,940 which is an ELMy H-mode case. Magnetic equilibrium and profiles are constructed prior to ELM onset. This particular case has been used recently for micro-tearing mode studies [27,28]. Details about this particular case can be found in [27], and further experimental details regarding DIII-D discharge 162,940 can be found in [29]. We use a Miller equilibrium [30] and obtain the Miller parameters from the EFIT equilibrium with a $513\times 513$ $(R,Z)$ grid, along with the density and temperature profiles.

The purpose of this work is to directly compare theoretical models, not to predict experimental transport levels. We do not include the effect of the equilibrium shear flow, as simple quasi-linear theory does not take this into account. Including the zonal flow and cross-coupling between the electron and ion scales continue to be active research topics [4,11]. Neglecting the shear flow in the quasilinear theory allows for a more transparent comparison of the various saturation rules. We include realistic collisionality in both the linear analysis and the nonlinear gyrokinetic simulations. Gyrokinetic ions and drift-kinetic electrons with electromagnetic fluctuations perpendicular to B ($\delta B_{\perp }$) are used in the linear and nonlinear gyrokinetic simulations, while $\delta B_{\Vert }$ is neglected. The plasma $\beta$ is reduced for a number of the nonlinear simulations presented in Section 5, the details of which are discussed there.

Figure 1a shows the main ion (Deuterium) and electron density profiles for DIII-D 162,940, Figure 1b shows the carbon impurity density profile, Figure 1c shows the electron temperature profile, and Figure 1d shows the main ion temperature profile. For this study, the impurity temperature is assumed to be equal to the main ion temperature. We only account for one impurity species, namely, carbon. The carbon profile is hollow; hence, we expect inward radial flux for the impurity species. We chose the three radial locations shown in

Table 1. Local tokamak plasma parameters at $\rho$ = 0.8, 0.85, and 0.9.

$\mathit{\rho }$	$\frac{\mathit{R}}{{\mathit{L}}_{{\mathit{T}}_{\mathit{i}}}}$	$\frac{\mathit{R}}{{\mathit{L}}_{{\mathit{T}}_{\mathit{e}}}}$	$\frac{{\mathit{T}}_{\mathit{e}}}{{\mathit{T}}_{\mathit{i}}}$	$\frac{\mathit{R}}{{\mathit{L}}_{\mathbf{ne}}}$	$\frac{\mathit{R}}{{\mathit{L}}_{\mathbf{ni}}}$	$\frac{\mathit{R}}{{\mathit{L}}_{\mathbf{nC}}}$	$\frac{{\mathit{n}}_{\mathit{C}}}{{\mathit{n}}_{\mathit{e}}}\mathbf{[}\mathbf{\%}\mathbf{]}$	$\mathit{q}$	$\hat{\mathit{s}}$	${\mathit{\beta }}_{\mathit{e}}\mathbf{[}\mathbf{\%}\mathbf{]}$	$\mathit{\kappa }$	$\mathit{\delta }$	$\mathit{\zeta }$
0.8	6.71	7.49	0.87	1.44	2.40	−0.71	5.16	2.28	1.75	1.00	1.47	0.21	−0.03
0.85	9.16	8.94	0.87	1.84	3.06	−0.75	5.37	2.56	2.17	0.85	1.51	0.24	−0.04
0.9	12.49	11.17	0.88	2.36	3.98	−0.79	5.65	2.97	2.94	0.67	1.55	0.28	−0.05

for our study, where $\rho =r/a$ and r is the Miller radial coordinate [30]:

$$r=\frac{R_{\max }-R_{\min }}{2},$$

(1)

where $R_{\max }$ and $R_{\min }$ are the maximum and minimum major radii of each flux surface, respectively, and a is the value of r from Equation (1) at the separatrix.

Table 1 provides the local parameters at the three radial locations ($\rho$) used in our analysis. We begin by examining the linear stability at these three radial locations. Note that the quasilinear analysis in Section 4 is a local analysis and is based on the local parameters provided in Table 1. Physical quantities such as the major radius, B, n, and T are important for determining collisionality and for conversion to physical (SI) units. Here, $\frac{R}{L_{Ti}}$ and $\frac{R}{L_{Te}}$ are the normalized temperature gradients of ions and electrons, respectively. We assume the carbon temperature (and temperature profile) to be the same as the main ions, with $\frac{R}{L_{ni}}$, $\frac{R}{L_{ne}}$, and $\frac{R}{L_{nc}}$ being the normalized density gradients of ions, electrons, and carbon, respectively. The ratio of the electron temperature to the ion temperature and the impurity concentration (carbon density to electron density) are respectively provided by $\frac{T_e}{T_i}$ and $\frac{n_C}{n_e}$. Here, q is the safety factor, $\hat{s}=\frac{\rho }{q}\frac{dq}{d\rho }$ is the magnetic shear parameter, and ${\beta }_e={\mu }_0{n}_e{T}_i/B^2$, while The Miller parameters [30] for elongation, triangularity, and squareness are $\kappa$, $\delta$, and $\zeta$, respectively.

3. Linear Analysis

We begin by studying the local linear properties of the tokamak plasma parameters (162,940) discussed above in Section 2 near the pedestal top, and scan $\rho =0.8,0.85,0.9$. We perform initial value calculations with the GENE code in the flux tube limit. Figure 2 shows the linear growth rate and real frequency versus $k_y$ for the three radial locations specified in Table 1, where y is the binormal perpendicular coordinate. In the following section, we use the linear output from GENE in the form of the particle and energy fluxes, along with the electrostatic potential linear mode structure, to parameterize $k_{\perp }$.

Here, ${\rho }_i$ is the deuterium species ion gyroradius and $v_{th}=\sqrt{T_i/m_i}$. Throughout this paper, the subscript “i” refers to the main deuterium ion species. Figure 2 shows the qualitative features common to core H-mode plasmas, and even the so-called “cyclone base case” [31]. An ion mode or ion temperature gradient (ITG) mode dominates for $k_y{\rho }_i\lesssim 1.4$, while an electron mode or collisionless trapped-electron mode (CTEM) dominates for $k_y{\rho }_i\gtrsim 1.4$. At $\rho =$ 0.8 and 0.85, there is an unstable mode with a negative frequency at the longest resolved wavelength. This is the micro-tearing mode (MTM), which is often the dominant mode in the pedestal region for these parameters [27,32]. Global analysis is required to accurately model the long wavelength (MTM).

Figure 3 shows a linear comparison between TGLF and local GENE at $\rho =$ 0.85. In Section 4, we compare the quasilinear theory results to TGLF as well. The real frequencies between GENE and TGLF agree very well. Note that the TGLF model [4] predicts higher growth rates than GENE. While GENE provides more accurate local linear gyrokinetic calculation, it is possible for the higher growth rate to be compensated by the TGLF saturation rule, i which case TGLF provides accurate fluxes. Note that global effects are typically be stabilizing, making this another important effect for realistic modeling that is not accounted for here.

4. Quasilinear Theory

4.1. Quasilinear Expression of Fluxes Using Linear Gyrokinetic Simulation

In good confinement regimes, core tokamak turbulence fluctuations are small. It is not unreasonable to assume a superposition of a finite number of linear eigenmodes at small amplitude, leading to the validity of the quasilinear expression for the fluxes. The quasilinear flux is quadratic in the mode amplitude. What is more uncertain (or unknown) is the fluctuation amplitude; we discuss three plausible saturation rules in the following section. Linear flux tube gyrokinetic simulation is used to predict the quasilinear fluxes by assuming a saturation rule. We follow Lapillonne’s prescription [15] and use the GENE code to obtain the linear fluxes. GENE uses field line-following coordinates $(x,y,z)$, where x is a radial coordinate, y is the other binormal coordinate, and z is the coordinate along the field line. We decompose the fluxes in $k_y$ and define a general quasilinear flux quantity $F^{ql}$, where F can represent the particle flux ${\Gamma }_{\alpha }$ or energy flux $Q_{\alpha }$ for species $\alpha$. The linear flux is proportional to the square of the amplitude

$$F_{k_y}^{lin}={\hat{G}}_{k_y}{\left|{\hat{\mathsf{\Phi }}}_{0,k_y}(z=0)\right|}^2,$$

(2)

where we assume that the mode amplitude can be parameterized at $z=0$ (or at the poloidal angle $\theta =0$). Equation (2) defines ${\hat{G}}_{k_y}$. It is straightforward to calculate the amplitude normalized linear flux ${\hat{G}}_{k_y}$ from linear gyrokinetic simulations. With a saturation rule for the amplitude, the fluxes can be calculated using

$$F^{ql}=\sum _{k_y}{A}^2\left(k_y\right){\hat{G}}_{k_y}^{ql}\Delta k_y,$$

(3)

where $\Delta k_y$ is the $k_y$ spacing while $A^2\left(k_y\right)$ parameterizes the mode amplitude and is determined using the three simple saturation rules discussed below. The saturation rule used in Lapillonne [15] is

$$A^2\left(k_y\right)=A_0^2{\left(\frac{{\gamma }_{k_y}}{〈k_{\perp }^2〉}\right)}^2,$$

(4)

where ${\gamma }_{k_y}$ is the linear growth rate. Equation (4) is one of the saturation rules we examine. Care must be taken in determining $〈k_{\perp }^2〉$ in the denominator of Equation (4), and we follow a similar procedure here to ensure consistency with previous work [15,21]. We average $k_{\perp }$ over the eigenmode envelope ${\hat{\mathsf{\Phi }}}_{k_x,k_y}\left(z\right)$ provided by

$$〈k_{\perp }^2〉=\frac{{\sum }_{k_x}\int (g^{xx}{k}_x^2+2g^{xy}{k}_x{k}_y+g^{yy}{k}_y^2){\left|{\hat{\mathsf{\Phi }}}_{k_x{k}_y}\left(z\right)\right|}^{2}Jdz}{{\sum }_{k_x}\int {\left|{\hat{\mathsf{\Phi }}}_{k_x{k}_y}\left(z\right)\right|}^{2}Jdz},$$

(5)

where J is the Jacobian and $g^{xx},g^{xy},g^{yy}$ are the geometric coefficients $g^{\mu \nu }=\nabla \mu ·\nabla \nu$ in the field line-following coordinates [15]. Note that $〈k_{\perp }^2〉$ in Equation (5) is a function of $k_y$.

4.2. Saturation Rules

Admittedly, although there may be validity in the quasilinear expression for the fluxes, the parameter dependence of the fluctuation amplitude is difficult to determine without running a nonlinear gyrokinetic simulation many times over a range of parameters [3]. However, we can gain some knowledge of the transport properties by comparing different saturation rules and the sensitivity of the observed trends. In this paper, we compare three common saturation rules and provide simple scaling arguments for their origin where they exist. The first saturation rule, in Equation (4), can be obtained by balancing the linear growth with the $\mathbf{E}\times \mathbf{B}$ advection. For example, the mode becomes saturated when $\frac{\partial \delta n}{\partial t}$ balances ${\mathbf{v}}_E·\nabla \delta n$, where ${\mathbf{v}}_E$ is the $\mathbf{E}\times \mathbf{B}$ drift, resulting in

$$\gamma \delta n_k\sim \frac{k_{\perp }^2}{B}\left|{\varphi }_k\right|\left|\delta n_k\right|$$

or

$$\frac{e\left|{\varphi }_k\right|}{T}\sim \frac{eB}{T}\frac{{\gamma }_k}{k_{\perp }^2},$$

which is the scaling in Equation (4), and is commonly used [14,15,16]. The saturation rule can be obtained from wave particle trapping of resonant particles, and provides the correct saturation level in slab geometry [33].

The second saturation rule we use comes from a dimensional argument, where the diffusion coefficient is simply set to $D=\gamma /k_{\perp }^2$. Then, $D\nabla n_0=〈v_{Ex}\delta n〉$, and we obtain

$$\frac{e{\left|{\varphi }_k\right|}^2}{T}=A_0^2\frac{eB}{T}\frac{1}{L{k}_y}\frac{{\gamma }_k}{k_{\perp }^2},$$

(6)

where L is the gradient scale length. In this argument, $L=L_n$, the density gradient scale length. A similar calculation can be made for the thermal diffusivity; thus, we write L in Equation (6) more generally. Equation (6) is similar to the saturation rule used in QuaLiKiz [9] as well as in earlier work [13].

Finally, the third saturation rule we use for comparison is

$$\frac{e{\left|{\varphi }_k\right|}^2}{T}=A_0^2\frac{eB}{T}\frac{{\gamma }_k}{k_{\perp }^2},$$

(7)

which has a similar $\frac{\gamma }{k_{\perp }^2}$ scaling as Equation (6) but does not diverge as $k_y$ approaches zero. The saturation rule in Equation (7) has been used previously for comparison to nonlinear gyrokinetic simulation and experiment [21]. Our goal here is not to develop a transport model that accurately reproduces nonlinear gyrokinetic simulation; rather, as we admit that the saturation rule is the weakness in any weak turbulence model, we present results for the three saturation rules above and in some sense “scan” the sensitivity of the fluxes to the saturation rule.

4.3. Saturation Levels and Quasilinear Fluxes

Figure 4 shows the three saturation rules in Equations (4), (6) and (7) with $A_0=1$ and $T=T_i$ along with the TGLF SAT1 results. GENE gyroBohm units are used, in which $\varphi$ is normalized by $\frac{R}{{\rho }_i}\frac{e\varphi }{T_i}$, i.e., to convert to SI units [V²], we take the values presented in Figure 4 and multiply by ${\left(\frac{{\rho }_i}{R}\frac{T_i}{e}\right)}^2$. The saturation rules in Equations (4), (6) and (7) are labeled “Lapillonne(2011)”, “Bourdelle(2007)”, and “Kumar(2021)”, respectively, for the convenience of the reader. TGLF SAT1 is labeled “SAT1”. No indication of the relative validity of the various models should be made here, as we are only using the various saturation rules for comparison. Additionally, the theory-based transport models calculate the fluxes differently than we do here, as we use the linear GENE results. The overall level of each saturation rule, i.e., the value of $A_0$, has little meaning, as quasilinear transport models calibrate the saturation rule using a constant coefficient.

The saturation rules show similar trends, peaking at $k_y{\rho }_i\sim$ 0.2–0.3, with variation in the width of the spectra. SAT1 has the narrowest spectrum, and not surprisingly Equation (7) (Kumar(2021)) has a broader spectrum compared to Equation (4) (Lapillonne(2011)) and Equation (6) (Bourdelle(2007)) due to the power of the $\frac{{\gamma }_k}{k_{\perp }^2}$ term and the lack of a $\frac{1}{k_y}$ factor. It is interesting that Equation (6) and SAT1 are somewhat similar.

Next, we compare the quasilinear flux obtained from the normalized GENE quasilinear flux and the three saturation rules. The results are shown in Figure 5. We use GENE gyroBohm units, where $Q\left[\mathrm{SI}\right]=\left(\frac{v_{th}{\rho }_i^2{n}_e{T}_i}{R^2}\right)Q_{\mathrm{shown}}$ and $\Gamma \left[\mathrm{SI}\right]=\left(\frac{v_{th}{\rho }_i^2{n}_e}{R^2}\right){\Gamma }_{\mathrm{shown}}$. The results in Figure 5 have been normalized by adjusting $A_0$ in the three saturation rules to ensure that the total ion heat flux $Q_i$ matches that predicted by SAT1. The SAT1 result shown in Figure 5 is from using the TGLF model directly.

The carbon flux is directed inward, as expected, due to the slightly hollow profile of carbon. The results from all four models agree qualitatively, again with some variation in the breadth of the spectrum, with TGLF showing more flux at lower $k_y$. The relatively large values of TGLF fluxes versus $k_y$ is simply an artifact of normalizing the other quasilinear fluxes to TGLF and the fact the other models have a broader $k_y$ spectrum. TGLF shows some electron flux at higher $k_y$ that the GENE quasilinear fluxes do not.

5. Comparison with Nonlinear Gyrokinetic Simulations

Local nonlinear flux tube simulations were carried out using the GENE code at $\rho =0.85$ to test the validity of the quasilinear models. The perpendicular flux tube domain we used was 167${\rho }_i$ (radial) × 126${\rho }_i$ in size, with 256 radial grid points and 32 toroidal modes and with $k_y{\rho }_i$ ranging from $0.05$ to $1.60$. The grid resolution in the $z,v_{\Vert },$ and $\mu$ dimensions was chosen to be $32\times 32\times 16$, respectively, with the values found by running linear growth rate convergence tests. Initial runs with the value of ${\beta }_e=0.85\%$ taken from Table 1 showed nonlinearly excited lower-$k_y$ micro-tearing modes (MTM) dominating the transport at earlier times, e.g., $t\frac{v_{th}}{R}\le 20$, and eventually becoming numerically unstable at later times. Because electromagnetic modes were not observed in the global nonlinear GEM simulations discussed below, ${\beta }_e$ was reduced to 10% of the original value, with high-quality electrostatic simulation results obtained in consequence. This change is reasonable as the fluxes are mainly electrostatic in the GEM simulations, which is because nonlocal effects can help to stabilize the low-n electromagnetic modes.

Figure 6 shows a comparison of the particle and heat fluxes versus $k_y$ between local nonlinear GENE, quasilinear theory, and global nonlinear GEM. The quasilinear fluxes shown here use linear GENE in combination with the saturation rule provided in Equation (4). The quasilinear fluxes shown in blue and labeled “QL GENE” are in good agreement with the nonlinear GENE results, shown in orange and labeled “NL flux-tube GENE.” The amplitude $A_0$ is normalized using the nonlinear GENE ion heat flux. This is appropriate, as the value of $A_0$ is undetermined in the theory. The global GEM results discussed below are scaled by a factor of $3.49$. Nonlocal effects, including the profile, q, and magnetic shear variation are stabilizing. Therefore, it is typical that global calculations are more stable, producing lower fluxes.

The results labeled “Nonlinear global GEM” in Figure 6 are nonlocal and nonlinear electromagnetic gyrokinetic simulations using the $\delta f$ particle-in-cell code GEM. For the present study of ion-scale turbulent transport, a fully drift-kinetic electron species is included using the $\delta f$ method with an expansion around the electron adiabatic response [26]. To ensure steady-state turbulence and transport, a fixed-gradient heat source is applied to all species [34]. The grid resolution is $(N_x,N_y,N_z)=(128,128,64)$ in the radial, binormal, and parallel directions, respectively. The particle number is 32/cell for the ion species and 64/cell for electrons. The time step is ${\mathsf{\Omega }}_p▵t=1$, where ${\mathsf{\Omega }}_p$ is the proton gyrofrequency. The radial domain of the nonlocal simulation is $0.65<r/a<0.95$. Attempts to extend the simulation to the separatrix ($r/a=1.0$) lead to nonphysical modes near the edge. While the cause of this problem is not clear, we believe part of the reason is the uncertainty in the equilibrium configuration, including the magnetic surface shape and the density/temperature profiles. In particular, strong poloidal variation of the temperature is expected in a region of steep gradients; however, such variation is not modeled in the present study, as local Maxwellian distributions that vary only in the radial flux coordinate are assumed. The density and temperature gradients are reduced in a boundary layer of $▵r/a\sim 0.05$ near the outer domain to avoid peaking of the turbulence near the boundary.

In the toroidal direction, the simulation domain is a toroidal wedge which is $1/8$ of the torus, and the EM fields are filtered to include only the toroidal mode numbers $n=0,8,16,\dots ,88$. We note that no nonlinear excitation of low-n electromagnetic modes, e.g., MTMs, was present in the global GEM simulations, possibly due to nonlocal stabilization, in contrast to local GENE results. Figure 7 shows the ion energy fluxes at various radial locations from GEM. The quality of the simulation seems reasonable.

For comparison to GENE and quasilinear theory, the turbulent fluxes are decomposed into toroidal modes. For example, in the formula for the radial heat flux

$$Q\left(r\right)=\frac{1}{▵V}{\int }_{▵V}\frac{1}{2}m{v}^2\delta f\left(\frac{\mathbf{E}\times \mathbf{B}}{B^2}+v_{\Vert }\frac{\delta {\mathbf{B}}_{\perp }}{B}\right)·\frac{\nabla r}{\left|\nabla r\right|}d\mathbf{x}d\mathbf{v},$$

if $\mathbf{E}$ and $\mathbf{B}$ are replaced by a specific toroidal component, then the contribution of that component to the total flux is obtained. Here, $▵V$ is a thin toroidal annulus with a radial size of $▵r/a=0.025$. Results from this procedure are shown in Figure 8. The raw results are shown as solid red triangles, while the solid blue squares are a polynomial fit. The fitted result is shown in Figure 7. GEM is a particle code, and does not evolve the distribution function spectrally in $k_y$ as the GENE calculation does. Obtaining $Q\left(k_y\right)$ involves summing over the particle weights for each toroidal mode in GEM, and leads to statistical fluctuations in this quantity. The smooth fit shown in Figure 6 (the solid blue squares in Figure 8) more clearly compares the trends between models. We carried out the same fitting procedure for all the flux quantities in Figure 6.

6. Discussion and Future Work

In high-confinement tokamak regimes, it is reasonable to assume that the quasilinear expression for the flux is valid; however, there is uncertainty in determining the fluctuation level. Here, we have compared three common saturation rules using gyrokinetic simulation to determine the quasilinear flux expression. To the best of our knowledge, the various common saturation rules have not been compared in this way before. Considerable realism was taken into account, including experimental plasma parameters, collisionality, electromagnetics, and drift-kinetic electrons. Quasilinear theory in combination with a chosen saturation rule can only predict parametric dependence of the transport level. Thus, an arbitrary constant (or possibly a slowly varying function), e.g., $A_0$ in Equation (4), needs to be determined by other means. Here, we show that with proper normalization to nonlinear gyrokinetic simulation or experiment, the three saturation rules exhibit similar behavior. We additionally compared local and global gyrokinetic simulations. The local nonlinear GENE results exhibited nonlinearly excited low-$k_y$ electromagnetic modes that were dominant. However, when $\beta$ was reduced in the simulation to eliminate the low-$k_y$ modes, nonlinear GENE showed similar behavior to quasilinear theory. Moreover, electromagnetic global GEM showed qualitatively similar behavior (see Figure 6). GEM did not see a low-n electromagnetic mode, possibly due to nonlocal stabilization, e.g., variation in profiles, q, and magnetic shear in this region. Global GEM provided a lower ion and electron particle flux, while nonlinear GENE and GEM showed higher impurity fluxes. All calculations presented in this paper are quasineutral, of course. Our further work will include a better understanding of the larger inward particle flux found in the nonlinear simulations.

This study has only scraped the tip of the iceberg, and further investigation is needed in order to better understand the parametric dependence of the models for a wider variety of plasma equilibria and profiles. Even for this particular case (DIII-D 162940), it would be interesting to look at the effect of including the equilibrium shear flow and the radial dependence of the fluxes. These are topics for future work. Detailed comparisons with nonlinear simulations are challenging due to computing resource restrictions. However, comparison between the various quasilinear models is relatively fast, requiring only linear GENE calculations. Additionally, it would be useful to investigate the relative contributions of particle diffusion and of the convection and impurity peaking factor [35,36]. Because the amplitude cancels out in the determination of the peaking factor in quasilinear theory, the choice of the particular saturation rule should be less important.