2.1. Linear simulation results

Faithfully capturing the linear physics will be critical to the quasi-linear model, so here we briefly revisit the linear analysis of STEP-EC-HD, detailed in Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023). Linear local gyrokinetic simulations were carried out on six radial surfaces using the GS2 code (Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Barnes et al. Reference Barnes, Dickinson, Dorland, Hill, Parker, Roach, Biggs-Fox, Christen, Numata and Jason2022). These calculations included three plasma species (electron, deuterium and tritium) and used the following numerical resolutions: extended parallel domain in ballooning angle spanning $-7{\rm \pi} <\theta < 7 {\rm \pi}$ ($n_\mathrm {period}=4$), 65 parallel grid points in each [$-{\rm \pi}, {\rm \pi}$] interval ($n_\theta =65$), 40 pitch angles ($n_\lambda = 40$) and 12 energy grid points ($n_\epsilon = 12$). The extended ballooning domain is increased to $n_\mathrm {period}=32$ on the innermost surface at $\varPsi _n=0.15$ in order to resolve more extended MTMs found on this surface. All calculations are fully electromagnetic, including $\delta \phi$, $\delta A_\parallel$ and $\delta B_\parallel$,Footnote ³ and were carried at $\theta _0=0$ unless otherwise specified.

Figure 1 shows the growth rate and the mode frequency as functions of the binormal wave vector $k_y\rho _s$ on the six surfaces. No unstable modes are found at $k_y\rho _s>1.8$ on any of the surfaces considered here. The maximum normalised local growth rate generally increases from the innermost to the outermost radial surface, and a sign change in the dominant mode frequency from mostly negative across most of the spectrum at $\varPsi _n=0.15$ to mostly positive at $\varPsi _n=0.36$, clearly indicates an instability transition going outwards from the core.

Figure 1.

Growth rate (a) and mode frequency (b) as functions of the binormal wave vector $k_y$ at various radial surfaces of STEP-EC-HD. The inset in (b) shows the transition to positive frequency values at low $k_y$ of $\varPsi _n = 0.15$.

The dominant mode at $\varPsi _n=0.49$ has been labelled as a hybrid-KBM because it shares many properties with KBM and exhibits ITG and TEM contributions to the linear drive (Kennedy et al. Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023). Figure 2 shows the parallel mode structure of $\delta \phi$ and $\delta A_\parallel$ of the dominant hybrid-KBM at $k_y\rho _s\simeq 0.4$ on $\varPsi _n=0.49$.

Figure 2.

Parallel mode structure of $\delta \phi$ (a) and $\delta A_\parallel$ (b) both normalised to $\max (\delta \phi )$ for the most unstable mode with $k_y\rho _s \simeq 0.4$. Figure adapted from Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023).

On the innermost surface at $\varPsi _n=0.15$, however, the dominant modes have negative frequency and very different parallel mode structures. These are illustrated in figures 3(a) and 3(b) for the fastest growing mode at $k_y\rho _s \simeq 0.84$: $\delta \phi$ has odd parity and is highly extended in $\theta$ whilst $\delta A_\parallel$ is very localised at $\theta =0$ with even parity. Amplitudes of $e\delta \phi /T_e$ and $\delta A_\parallel /(\rho _s B_0)$ are comparable, indicating that the mode is strongly electromagnetic. These negative frequency modes propagate in the electron diamagnetic drift direction, and have the typical properties expected of MTMs (Applegate et al. Reference Applegate, Roach, Cowley, Dorland, Joiner, Akers, Conway, Field, Patel and Valovic2004). On the same surface at the lowest $k_y$, on the other hand, the dominant mode switches to a positive frequency (see the inset in figure 1b) and very different parallel mode structures in $\delta \phi$ and $\delta A_\parallel$. The parallel mode structure of a low $k_y$ mode with positive frequency is shown in figures 3(c) and 3(d): the parallel structure closely resemble figure 2, indicating that the hybrid-KBM is linearly dominant at low $k_y$. Although these low $k_y$ hybrid-KBMs are only weakly unstable, they nevertheless contribute significantly to the turbulent transport, as will be discussed in § 2.2.

Figure 3.

Parallel mode structure of $\delta \phi$ (a,c) and $\delta A_\parallel$ (b,d) both normalised to $\max (\delta \phi )$ for the unstable mode with $k_y\rho _s \simeq 0.84$ (a,b) and $k_y\rho _s \simeq 0.05$ (c,d) at $\varPsi _n=0.15$.

The analysis of Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023) showed that hybrid-KBM is highly ballooned with growth rates that are extremely sensitive of the ballooning parameter, $\theta _0$. This suggests that hybrid-KBM turbulence should be sensitive to equilibrium flow shear, as was subsequently confirmed in Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024). Figure 4(b) shows the dependence on $\theta _0$ of growth rates on the outermost surface at $\varPsi _n=0.8$, where hybrid-KBMs similar to the modes at $\varPsi _n=0.49$ dominate: the hybrid-KBM growth rates are strongly suppressed as $\theta _0$ increases from zero (the slowly growing unstable mode on this surface at $\theta _0={\rm \pi}$ and $k_y\rho _s \simeq 0.03$ is a MTM). Figure 4(a), however, shows that the situation is quite different on the innermost surface at $\varPsi _n=0.15$: here, the MTM growth rates at $k_y\rho _s \simeq 0.8$ and $k_y\rho _s \simeq 1.3$ are insensitive to $\theta _0$,Footnote ⁴ although in contrast the dominant hybrid-KBM at $k_y\rho _s\simeq 0.05$ has a growth rate that decreases with $\theta _0$, albeit remaining unstable across the entire $\theta _0$ domain.

Figure 4.

Dependence of the growth rate on $\theta _0$ at three different values of $k_y\rho _s$ on the surfaces at $\varPsi _n=0.15$ (a) and $\varPsi _n = 0.8$ (b). Stable modes are shown with $\gamma = 0$ and open markers.

2.2. Nonlinear simulations results

This section presents a brief overview of local nonlinear gyrokinetic simulations carried out on three radial surfaces, $\varPsi _n\in \{0.15,\;0.36,\; 0.8\}$, of STEP-EC-HD. This analysis extends the investigation carried out at $\varPsi _n=0.49$ by Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024), providing further nonlinear simulations to be used in the development of a reduced turbulent transport model to account for hybrid-KBMs. The nonlinear gyrokinetic simulations are carried out with the GENE code (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000) using the advanced Sugama collision operator (Sugama, Watanabe & Nunami Reference Sugama, Watanabe and Nunami2009). The numerical resolutions used in these simulations is reported in table 2. We highlight the higher radial resolution required by MTMs at $\varPsi _n=0.15$. Simulations are performed without and with equilibrium flow shear. All the simulations are performed over a time interval that covers at least a growth time of the slowest growing low $k_y$ mode of each surface, corresponding to several growth times of the most unstable mode.

Table 2.

Numerical resolution of the nonlinear GENE simulations at the different radial locations. In the table, $n_\mu$ is number of grid points in the $\mu =v_{\perp }^2/(2B)$ direction and $n_v$ is number of velocity grid points.

2.2.1. Innermost surface dominated linearly by MTMs ($\varPsi _n=0.15$)

Figure 5 shows the time trace of the total heat and particle flux $k_y$ spectrum, normalised to $Q_\mathrm {gB} = \rho _*^2 n_e T_e c_s$ and $\varGamma _\mathrm {gB} = \rho _*^2 n_e c_s$, from a nonlinear simulation at $\varPsi _n=0.15$ without equilibrium flow shear. The heat flux contribution from $k_y\rho _s \sim 0.25$ increases until $t\simeq 250\,a/c_s$ where it reaches values of the order of $Q_\mathrm {gB}$, and then remains approximately constant until $t\simeq 600\,a/c_s$. The particle flux contributions from $k_y\rho _s > 0.2$ are negligible, as expected from MTM-driven turbulent transport in this $k_y$ range on this surface (see, e.g. Doerk et al. (Reference Doerk, Jenko, Pueschel and Hatch2011), Guttenfelder et al. (Reference Guttenfelder, Candy, Kaye, Nevins, Wang, Bell, Hammett, LeBlanc, Mikkelsen and Yuh2011), Hatch et al. (Reference Hatch, Kotschenreuther, Mahajan, Valanju, Jenko, Told, Görler and Saarelma2016) and Giacomin et al. (Reference Giacomin, Dickinson, Kennedy, Patel and Roach2023) for details on MTM turbulence). On the other hand, the heat and particle flux contributions from $k_y\rho _s < 0.2$ steadily grow (the total heat flux grows to exceed 100 $Q_\mathrm {gB}$) and show no signs of saturation within the simulated time interval. Figure 5 suggests that the large turbulent fluxes at low $k_y$ are driven by the hybrid-KBMs that are linearly unstable at the lowest binormal wavenumbers in figure 1.

Figure 5.

Total heat (a) and particle (b) flux as functions of $k_y$ and time from the nonlinear simulation at $\varPsi _n=0.15$. The equilibrium flow shear is not included.

The sensitivity of $\gamma ^{\rm KBM}(k_y \rho _s=0.05)$ to $\theta _0$ shown in figure 4(a), together with the low hybrid-KBM growth rate, implies that equilibrium flow shear may be important here even at small shearing rates near the diamagnetic level, $\gamma _E^{\rm dia}$, where

(2.1)\begin{equation} \gamma_E^\mathrm{dia} = \frac{\rho}{q}\frac{\mathrm{d}}{\mathrm{\rho}}\left(\frac{E_r}{R B_\theta}\right) = \frac{\rho}{q}\frac{\mathrm{d}}{\mathrm{\rho}}\left[\frac{1}{R B_\theta}\left(\frac{1}{n_i e}\frac{\mathrm{d}p_i}{\mathrm{d}r}+v_\phi B_\theta - v_\theta B_\phi\right)\right], \end{equation}

with $E_r$ the neoclassical radial electric field, $n_i$ and $p_i$ the thermal ion density and pressure, $v_\theta$ and $v_\phi$ the poloidal and toroidal velocities and $B_\theta$ and $B_\phi$ the poloidal and toroidal magnetic field components. In the JINTRAC flat-top operating point considered here, $v_\theta$ is evaluated using NCLASS (Houlberg et al. Reference Houlberg, Shaing, Hirshman and Zarnstorff1997) while $v_\phi =0$. Figure 6 shows the time trace of the electrostatic and electromagnetic heat and particle fluxes from a nonlinear simulation at $\varPsi _n=0.15$ with $\gamma _E = \gamma _E^\mathrm {dia} \simeq 0.01\, c_s/a$. When the equilibrium flow shear is included at $t=150\,a/c_s$, the heat flux reduces slightly and saturates at approximately $5\ {\rm MW}\ {\rm m}^{-2}$. The runaway fluxes are completely avoided in the simulation with $\gamma _E =0.01\, c_s/a$. Figure 7 shows, however, that, even with this level of flow shear, the saturated heat and particle fluxes nevertheless remain above the values required to achieve a transport steady state on this surface in STEP-EC-HD with the available heat and particle sources.

Figure 6.

Time trace of the electromagnetic and electrostatic heat (a) and particle (b) fluxes from the nonlinear simulation at $\varPsi _n=0.15$ with $\gamma _E=0.01\,c_s/a$. The equilibrium flow shear is active from $t=150\,a/c_s$, as indicated by the dashed vertical line.

Figure 7.

Total heat (a) and particle (b) fluxes from nonlinear simulations at $\varPsi _n=0.15$ with (red) and without (blue) equilibrium flow shear. The heat and particle fluxes in the simulation without equilibrium flow shear exceed the value of $10^2\ {\rm MW}\ {\rm m}^{-2}$ and $10^{22}\ \text {particles}\ ({\rm m}^2\ {\rm s})^{-1}$ with no saturation within the simulation time considered here, as indicated by the triangle on the top of the bar. The dashed horizontal line denotes the target heat and particle flux values at $\varPsi _n=0.15$ computed from the JETTO heat and particle sources.

2.2.2. Further out radially where hybrid-KBMs dominate ($\varPsi _n=0.36$)

Figure 8 compares the heat and particle fluxes averaged over time in the saturated phase from nonlinear simulations without and with equilibrium flow shear of the order of the diamagnetic level, $\gamma _E^\mathrm {dia}\approx 0.05\, c_s/a$. Although this nonlinear simulation at $\gamma _E=0$ appears to saturate within the time simulated, the corresponding heat and particle fluxes exceed STEP-EC-HD transport steady-state values for this surface by more than three orders of magnitude. When flow shear is included, both heat and particle fluxes reduce significantly, but remain well above the transport steady-state values, as shown in figure 8.

Figure 8.

Total heat (a) and particle (b) fluxes from nonlinear simulations at $\varPsi _n=0.36$ with (red) and without (blue) equilibrium flow shear. The dashed horizontal line denotes the target flux value at $\varPsi _n=0.36$ computed from JETTO heat and particle sources.

The outermost flux surface, where hybrid-KBMs also dominate ($\varPsi _n=0.8$):

Nonlinear simulations deliver similar findings at $\varPsi _n=0.8$. Figure 9 shows the heat and particle fluxes from nonlinear simulations without and with equilibrium flow shear ($\gamma _E=0.1\,c_s/a$). The simulation without equilibrium flow shear does not saturate within the considered simulation time interval, reaching values larger than $10^3\ {\rm MW}\ {\rm m}^{-2}$. The simulation with equilibrium flow shear saturates at much lower values, but turbulent fluxes remain an order of magnitude larger than the transport steady-state values for STEP-EC-HD at $\varPsi _n=0.8$.

Figure 9.

Saturated value of the total heat (a) and particle (b) fluxes from nonlinear simulations at $\varPsi _n=0.80$ with (red) and without (blue) equilibrium flow shear. The heat and particle fluxes in the simulation without equilibrium flow shear exceed $10^3\ {\rm MW}\ {\rm m}^{-2}$ and $10^{23}\ {\rm particles}\ ({\rm m}^2\ {\rm s})^{-1}$ with no saturation within the simulation time considered here, as indicated by the triangular arrow at the top of the bar. The dashed horizontal line denotes the target flux value at $\varPsi _n=0.80$ computed from the JETTO heat and particle sources.

3.1. Quasi-linear reduced transport model

In the quasi-linear inspired reduced model developed here, the saturation rule describes the total heat flux,Footnote ⁵ which is modelled asFootnote ⁶

(3.1)\begin{equation} Q_\mathrm{ql} = Q_0 \varLambda^\alpha, \end{equation}

where the quasi-linear metric $\varLambda$ (involving $\gamma /k_\perp ^2$, all fields and describing flow shear) is at the heart of the model, and $Q_0$ and $\alpha$ are constants chosen to fit the total turbulent heat transport in the regime of interest. In § 3.2, $Q_0$ and $\alpha$ will be optimised to describe the total heat flux from hybrid-KBM turbulence in STEP. Firstly we define $\varLambda$ and then we specify the species decomposition of the model transport fluxes.

In ballooning space local perturbed modes are labelled by the bi-normal perpendicular wavenumber and ballooning parameter, $k_y$ and $\theta _0$, respectively. For each mode (i.e. $k_y$ and $\theta _0$) and field $\chi _i \in \{e\delta \phi /T_e, \delta A_\parallel /(\rho _s B_0), \delta B_\parallel /B_0\}$, we define a perpendicular wavenumber averaged in poloidal angle $\theta$, following Jenko, Dannert & Angioni (Reference Jenko, Dannert and Angioni2005), as

(3.2)\begin{equation} \langle k_\perp^2\rangle_{k_y, \theta_0, i} = \frac{\displaystyle\int \mathrm{d}\theta\,k_\perp^2(\theta, k_y, \theta_0) \,J(\theta) \,|\chi_i(\theta)|^2}{\displaystyle\int \mathrm{d}\theta\, J(\theta)|\chi_i(\theta)|^2}, \end{equation}

where $k_\perp ^2(\theta, k_y, \theta _0) = g^{yy}(\theta )k_y^2 + 2g^{xy}(\theta ) k_x k_y + g^{xx}(\theta )k_x^2$, $k_x = \hat {s} \theta _0 k_y$, $g^{yy}$, $g^{xx}$ and $g^{xy}$ are metric tensor elements and $J(\theta )$ is the Jacobian (these geometric quantities are defined for GS2 in Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000).

By using $\langle k_\perp ^2\rangle _{k_y, \theta _0, i}$ from (3.2) and summing over all three fields, we define

(3.3)\begin{equation} \hat{\varLambda}(k_y, \theta_0) = \sum_{i=1}^3 \frac{\max_\theta|\chi_i(k_y,\theta_0, \theta)|}{\max_\theta| \chi_1(k_y,\theta_0, \theta) |} \,\frac{\gamma(k_y, \theta_0)}{\langle k_\perp^2(k_y, \theta_0, \theta)\rangle_{\theta,i}}. \end{equation}

When modes are strongly electromagnetic, $\hat {\varLambda }(k_y, \theta _0)$ has substantial contributions from magnetic perturbations.Footnote ⁷ On the other hand, when $\delta A_\parallel$ and $\delta B_\parallel$ vanish (e.g. at low $\beta$), (3.3) reduces to the electrostatic model of Jenko et al. (Reference Jenko, Dannert and Angioni2005).

The dependence of $\hat {\varLambda }$ on $\theta _0$ in (3.3) can be exploited to capture the effect of equilibrium flow shearFootnote ⁸ on the quasi-linear turbulent transport. Equilibrium flow shear advects local modes in the ballooning parameter at a rate $\mathrm {d}\theta _0/\mathrm {d}t=\gamma _E/\hat {s}$, which is accounted for by averaging $\hat {\varLambda }(k_y, \theta _0)$ over a suitable $\theta _0$ interval

(3.4)\begin{equation} \bar{\varLambda}(k_y) = \frac{1}{\theta_{0,\mathrm{max}}(k_y, \gamma_E)}\int_0^{\theta_{0,\mathrm{max}}(k_y,\gamma_E)}\hat{\varLambda}(k_y, \theta_0)\, \mathrm{d}\theta_0 , \end{equation}

where

(3.5)\begin{equation} \theta_{0,\mathrm{max}} = \min\left(\frac{\gamma_E}{\hat{s}\gamma}, {\rm \pi}\right). \end{equation}

The $\theta _0$ interval is the advection range over one growth time, from the typically most unstable mode at $\theta _0=0$ to $\theta _{0,\mathrm {max}}=\gamma _E/(\hat {s} \gamma )$, with $\gamma$ the growth rate at $\theta _0=0$. In the low flow shear limit, $\gamma _E \ll \hat {s} \gamma$, $\bar {\varLambda }(k_y)$ reduces to $\hat {\varLambda }(k_y,0)$; in the high flow shear limit where $\gamma _E > \hat {s} \gamma {\rm \pi}$, the average is calculated over a reduced $\theta _0$ range $[ 0, {\rm \pi}]$, justified by periodicity. The $\theta _0$ average in (3.4) allows the quasi-linear model to capture flow shear suppression of turbulence at finite $\gamma _E$. We note that, if $\gamma (\theta _0)$ is narrowly peaked close to $\theta _0=0$, $\bar {\varLambda }(k_y)$ can be sharply reduced even at small values of $\gamma _E$. As shown in figure 4 and in Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023), the growth rate of hybrid-KBMs can peak strongly at $\theta _0=0$ in STEP.

The quasi-linear metric $\varLambda$ in (3.1) is then obtained by integrating $\bar {\varLambda }(k_y)$ in (3.4) over $k_y$, and normalised to $\rho _* c_s$

(3.6)\begin{align} \varLambda & = \frac{1}{\rho_* c_s}\int \mathrm{d}k_y \,\bar{\varLambda}(k_y) \nonumber\\ & =\frac{1}{\rho_* c_s}\int \mathrm{d}k_y\frac{1}{\theta_{0,\mathrm{max}}}\int_0^{\theta_{0,\mathrm{max}}}\mathrm{d}\theta_0\,\sum_{i=1}^3 \frac{\max_\theta|\chi_i(k_y,\theta_0, \theta)|}{\max_\theta|\chi_1(k_y,\theta_0, \theta)|} \frac{\gamma(k_y, \theta_0)}{\langle k_\perp^2(k_y, \theta_0, \theta)\rangle_{\theta,i}} . \end{align}

This quantity depends only on linear physics, so its evaluation requires only linear gyrokinetic simulations over the unstable spectrum in $(k_y, \theta _0)$. We highlight that $\varLambda$ has been derived rather generally without any specific assumption on the turbulence regime. The parameters $Q_0$ and $\alpha$, which appear in (3.1), are the only free parameters to influence the turbulence saturation level in this reduced model.

Species heat and particle fluxes are obtained (from (3.7) and (3.8) below) using the total model heat flux, $Q_\mathrm {ql}=Q_0\varLambda ^\alpha$ defined in (3.1) that sets the turbulence saturation amplitude, together with the quasi-linear weights $Q_{{l},s}/Q_{l}$ and $\varGamma _{{l},s}/Q_{l}$, where $Q_{{l},s}$ and $\varGamma _{{l},s}$ are the linear heat and particle fluxes for species $s$, and $Q_{l} = \sum _s Q_{{l},s}$ is the total linear heat flux; all linear fluxes are taken from the final time point in a converged linear simulation for a given mode (i.e. $k_y, \theta _0$). The quasi-linear weights apportion contributions to transport fluxes for each species for a given mode (depending on $k_y$ and $\theta _0$) must be included in the integrands of (3.7) and (3.8) to obtain the total species heat and particle fluxes

(3.7)\begin{gather} Q_{\mathrm{ql},s} = Q_0 \varLambda^{\alpha-1}\left(\frac{1}{\rho_* c_s}\right) \int\mathrm{d}k_y\,\frac{1}{\theta_{0,\mathrm{max}}}\int_0^{\theta_{0,\mathrm{max}}} \mathrm{d}\theta_0 \,\frac{Q_{{l},s}(k_y,\theta_0)}{Q_{l}(k_y,\theta_0)} \hat{\varLambda}(k_y, \theta_0) \end{gather}

(3.8)\begin{gather}\varGamma_{\mathrm{ql},s} = Q_0 \varLambda^{\alpha-1}\left(\frac{1}{\rho_* c_s}\right) \int\mathrm{d}k_y\,\frac{1}{\theta_{0,\mathrm{max}}}\int_0^{\theta_{0,\mathrm{max}}} \mathrm{d}\theta_0 \,\frac{\varGamma_{\mathrm{l},s}(k_y,\theta_0)}{Q_{\mathrm{l}}(k_y,\theta_0)} \hat{\varLambda}(k_y, \theta_0), \end{gather}

where we note that $\sum _s Q_{\mathrm {ql},s} = Q_\mathrm {ql}=Q_0\varLambda ^\alpha$.

3.2. Application to STEP-relevant regimes

The parameters $Q_0$ and $\alpha$, which link $\varLambda$ to the quasi-linear heat flux $Q_\mathrm {ql}$ via (3.1), are tuned here to optimise the description of the total heat transport from hybrid-KBM turbulence in STEP. Suitable values are determined by fitting the nonlinear saturated heat fluxes from a set of nonlinear gyrokinetic simulations for STEP-relevant local equilibria. This database contains nonlinear simulations performed using local parameters from various reference surfaces in STEP-EC-HD. The parameter scans reported in Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023) and Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024) already cover many of the most important local parameters, including pressure gradient, $\beta$, flow shear, $q$ and $\hat {s}$, to which hybrid-KBM turbulence is sensitive. Those scans provide the bulk of the database used to develop the quasi-linear model. In detail, we include several parameter scans at $\varPsi _n=0.49$ with $\gamma _E\in \{0.05, 0.1, 0.2\}\,c_s/a$, $q\in \{3, 3.5, 4.0, 4.5\}$, $\beta \in \{0.005, 0.02, 0.09\}$ and $L_{p,\mathrm {ref}}/L_p\in \{0.8, 1.0, 1.2\}$ from Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024) and an additional scan with $\hat {s} \in \{0.3, 0.6, 1.2\}$ for $\gamma _E=0$. The database also includes the simulations at $\varPsi _n=0.36$ ($\gamma _E=0$ and $\gamma _E=0.05\,c_s/a$) and at $\varPsi _n=0.8$ ($\gamma _E=0.1\,c_s/a$) described in § 2.2, as well as an additional scan with $q\in \{3.0, 3.5, 4.0, 4.5\}$ for $\gamma _E=0$ at $\varPsi _n=0.36$.

Figure 10 compares the total normalised heat flux, $Q/Q_\mathrm {gB}$, from nonlinear simulations with that from the best-fit quasi-linear model defined in (3.1), as functions of $\varLambda$, where $\varLambda$ is defined in (3.6) and obtained from linear gyrokinetic simulations. Each point in figure 10 corresponds to a nonlinear simulation and the colour distinguishes simulations with and without equilibrium flow shear. The uncertainty in the average nonlinear heat flux from each simulation is taken as the standard deviation in the saturated phase.Footnote ⁹ These uncertainties are accounted for by using a weighted fit to obtain the best values for the model parameters $Q_0$ and $\alpha$. This gives $Q_0/Q_\mathrm {gB}\simeq 25$ and $\alpha \simeq 2.5$, with a fitting quality of $R^2\simeq 0.8$. The unweighted fit returns values of $\alpha$ and $Q_0$ that differ approximately by 20 %, which can be considered as an estimate of the uncertainty in the fit values. The quasi-linear model reproduces acceptably well the nonlinear heat flux across a range of more than three orders of magnitude, although we highlight significant scatter. The discrepancy between the reduced model and simulations approaches an order of magnitude in a few cases, illustrating a limitation in the reduced transport model. A comparison of the nonlinear species heat and particle fluxes with the quasi-linear predictions given by (3.7) and (3.8) is reported in the Appendix.

Figure 10.

Normalised total heat flux $Q_\mathrm {tot}/Q_\mathrm {gB}$ from nonlinear gyrokinetic simulations as a function of $\varLambda$. The dashed line represents the best weighted fit, $Q_\mathrm {tot}=Q_0 \varLambda ^\alpha$, with $Q_0$ and $\alpha$ the fitting coefficients. The error bar in the heat flux is taken as the standard deviation of the heat flux time trace in the saturated phase. Black and red markers correspond to simulations without and with equilibrium flow shear, respectively.

Turbulent transport in the nonlinear simulations of figure 10 is predominantly dominated by hybrid-KBM turbulence, but in two low $\beta$ cases turbulence is from electrostatic ITG and TEM. In principle, MTMs may generate magnetic stochastic transport if the separation between the magnetic islands from the tearing perturbation is smaller than the island width (see, e.g. Doerk et al. Reference Doerk, Jenko, Pueschel and Hatch2011; Guttenfelder et al. Reference Guttenfelder, Candy, Kaye, Nevins, Wang, Bell, Hammett, LeBlanc, Mikkelsen and Yuh2011; Nevins, Wang & Candy Reference Nevins, Wang and Candy2011; Giacomin et al. Reference Giacomin, Dickinson, Kennedy, Patel and Roach2023), but this physics is not captured here or in most other quasi-linear models since it would require predictions of the saturated magnetic island width generated by the underlying MTM instability. On the other hand, the nonlinear gyrokinetic simulations of Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024) suggest that magnetic stochastic transport should be modest on most STEP-EC-HD surfaces that have been analysed. Other simulations, however, suggest that stochastic transport may be significant in some spherical tokamak plasmas from NSTX (Guttenfelder et al. Reference Guttenfelder, Candy, Kaye, Nevins, Wang, Bell, Hammett, LeBlanc, Mikkelsen and Yuh2011) and MAST (Giacomin et al. Reference Giacomin, Dickinson, Kennedy, Patel and Roach2023). It is entirely possible that MTM induced magnetic stochastic transport will be substantial under some plasma conditions in STEP, and investigating this is an important area for future work. It is therefore likely that the reduced turbulent transport model developed here will need to be extended to account for magnetic stochastic turbulent transport from MTMs. We note that the development of such models is an active research area (Rafiq et al. Reference Rafiq, Weiland, Kritz, Luo and Pankin2016; Curie et al. Reference Curie2022; Hamed et al. Reference Hamed, Pueschel, Citrin, Muraglia, Garbet and Camenen2023; Hornsby et al. Reference Hornsby, Gray, Buchanan, Patel, Kennedy, Casson, Roach, Lykkegaard, Nguyen and Papadimas2024).

Figure 11(a) shows examples of the $\delta \phi$ and $\delta A_{\parallel }$ contributions to the $k_y$ spectrum of $\bar {\varLambda }(k_y)$ for the surface at $\varPsi _n=0.36$ of STEP-EC-HD without and with flow shear ($\gamma _E=0$ and $\gamma _E = 0.05\, c_s/a$). Both contributions to $\bar {\varLambda }(k_y)$ peak at low $k_y$ with and without flow shear, but at $\gamma _E=0$ the peaks are significantly larger and at longer wavelength. We note that, in this local equilibrium, the $\delta A_\parallel$ contribution to $\bar {\varLambda }(k_y)$ dominates both at $\gamma _E=0$ and $\gamma _E = 0.05\, c_s/a$, highlighting the need to include $\delta A_\parallel$ in (3.6).

Figure 11.

(a) Contributions of the $\delta \phi$ and $\delta A_\parallel$ terms in (3.4) to $\bar {\varLambda }(k_y)$ without (round markers) and with (squared markers) equilibrium flow shear. The contribution from $\delta B_\parallel$ is negligible and not shown. (b) Comparison between the nonlinear and the quasi-linear total heat flux for the two cases considered in (a). These simulations are performed at $\varPsi _n=0.36$ of STEP-EC-HD.

Figure 11(b) compares the nonlinear and quasi-linear model heat fluxes, where the reduced model used $Q_0$ and $\alpha$ in (3.1) from the fit in figure 11. The quasi-linear and nonlinear heat fluxes agree within 30 %, although it is pertinent to note that both simulations were included in the database used to obtain $Q_0$ and $\alpha$. The reduction in the heat flux attributable to flow shear approaches two orders of magnitude, and this is largely captured by the quasi-linear model. The reduction in the quasi-linear heat flux is slightly weaker than found in the nonlinear simulations, where the saturated total heat flux is 30 % lower at $\gamma _E=0.05\,c_s/a$.

We underline that the nonlinear heat flux calculations in the two cases of figure 11 requires more than $5 \times 10^5$ core hours for each simulation, corresponding to a week running on 32 nodes of the ARCHER2 high performance computing facility (Edinburgh, UK); in contrast computing the quasi-linear metric takes approximately 100 core hours, corresponding to 10 minutes running on 4 nodes of ARCHER2. Computing the quasi-linear metric is at least a factor of $10^3$ less expensive than computing the nonlinear heat flux.

Motivated by the analysis of Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024) that finds a significant influence of $q$ on hybrid-KBM turbulence at mid-radius in STEP-EC-HD, figure 12 explores the sensitivity of the quasi-linear heat flux to $q$ and $\hat {s}$ at $\varPsi _n=0.36$ and $\varPsi _n=0.49$. Equilibrium flow shear is not included in this two-dimensional parameter scan. At $\varPsi _n=0.49$, $Q_\mathrm {ql}/Q_\mathrm {gB}$ increases monotonically as $\hat {s}$ increases or $q$ decreases, which is consistent with the nonlinear $q$ scan presented in Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024). At $\varPsi _n=0.36$, the dependence of $Q_\mathrm {ql}/Q_\mathrm {gB}$ on $q$ and $\hat {s}$ is non-monotonic, with a maximum at $q=3.0$; at $q$ values higher or smaller than the nominal value, the quasi-linear heat flux decreases as $\hat {s}$ increases. In particular, figure 12 shows the presence of a low turbulent transport regime that could be achieved by varying $\hat {s}$ and $q$ from their nominal values, which are denoted by the white stars in figure 12; determining whether such local equilibrium conditions can be accessed across a significant radial region, while still satisfying operational constraints, clearly requires integrated scenario modelling (Tholerus et al. Reference Tholerus, Casson, Marsden, Wilson, Brunetti, Fox, Freethy, Hender, Henderson and Hudoba2024).

Figure 12.

Two-dimensional $(q, \hat {s})$ scan of the quasi-linear heat flux (normalised to $Q_\mathrm {gB}$) at $\varPsi _n=0.36$ (a) and $\varPsi _n=0.49$ (b). The white stars indicate the nominal values of $q$ and $\hat {s}$ of the chosen STEP flat-top operating point. The simulation at $\hat {s}=0.8$, $q=1.5$ and $\varPsi _n=0.36$ is stable. In this figure, the quasi-linear heat flux is computed without equilibrium flow shear.

4.1. Implementation of the reduced model in T3D

T3D (Qian et al. Reference Qian, Buck, Gaur, Mandell, Kim and Dorland2022) is a transport code that has been recently developed from TRINITY (Barnes et al. Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010). It is a very flexible and modular tool that facilitates the implementation of new transport models. The transport equations solved by T3D are (see Barnes et al. Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010)

(4.1)\begin{gather} \frac{\partial n_s}{\partial t} + \frac{1}{V'}\frac{\partial}{\partial \varPsi}(V'\overline{\varGamma_s}) = \overline{S_{n, s}}, \end{gather}

(4.2)\begin{gather}\frac{3}{2}\frac{\partial p_s}{\partial t} + \frac{1}{V'}\frac{\partial}{\partial \varPsi}(V'\overline{Q_s}) = \frac{3}{2}n_s\sum_u \nu_{su}(T_u - T_s) + \overline{S_{p,s}}, \end{gather}

where $n_s$ and $p_s$ are the density and pressure of the species $s$, $\varPsi$ is the flux label, $V$ is the volume enclosed by the surface $\varPsi$, $V'$ is the derivative of $V$ with respect to $\varPsi$, $\varGamma _s$ and $Q_s$ are the total particle and heat fluxes of a given species $s$ including the classical, neoclassical and turbulent contributions and $S_{n,s}$ and $S_{p,s}$ are the particle and the heat sources. The overline denotes a flux-surface average. The term proportional to $\nu _{su}$ is the collisional energy exchange term, with $\nu _{su}$ taken from Huba (Reference Huba1998). Density and temperature profiles are discretised on a uniformly spaced radial grid, and gradients are computed using a sixth-order finite difference scheme on interior radial grid points and a second-order finite difference scheme at the boundary. A non-zero Dirichlet boundary condition is applied at the outermost radial grid point, while the physical zero flux condition is applied at the magnetic axis.Footnote ¹⁰  Equations (4.1) and (4.2) are evolved using a semi-implicit method. Further details on the transport equations and on their numerical implementation are reported in Barnes et al. (Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010).

Evaluating the fluxes $\varGamma _s$ and $Q_s$ in (4.1) and (4.2) requires a transport model. For instance, the coupling between T3D and GX (Mandell et al. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2024), a local gyrokinetic code developed to run natively on graphics processing units, allows T3D to submit nonlinear gyrokinetic simulations at each radial location to compute $\varGamma _s$ and $Q_s$ from first-principles simulations. A total number of $(N_r-1)(N_p+1)$ simulations per iteration is required, where $N_r$ is the number of radial grid points and $N_p$ is the number of evolved profiles. We note that $N_p+1$ evaluations are required at each radial location to compute the local flux derivatives with respect to each evolved profile gradient (see Barnes et al. (Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010) for details). The quasi-neutrality condition is applied to impose the density value of the last evolved species. Our T3D simulations for STEP use the quasi-linear reduced transport flux equations (3.7) and (3.8), which have been implemented in T3D, and exploit a coupling to GS2 for the linear calculations required to compute the quasi-linear fluxes.

The workflow of T3D-GS2 is shown in figure 13. After initialisation, where the input file (containing run configuration data, such as the radial grid resolution, the maximum number of iteration and the transport model selection), geometry,Footnote ¹¹ initial profiles and sources are loaded, T3D-GS2 enters into the main loop, which is divided in four substeps. In the first substep, T3D writes the GS2 input files for each radial grid point and $(k_y, \theta _0)$ mode, according to the specified numerical resolutions. Several GS2 calculations are then distributed over the available MPI (Message Passing Interface) tasks. Depending on the numerical resolutions and on the number of linear gyrokinetic simulations to be performed, T3D-GS2 can efficiently scale to thousands of cores. Once all the GS2 linear simulations are completed, results are gathered and the quasi-linear fluxes are computed using (3.7) and (3.8). In the last sub-step, density and pressure profiles are evolved using (4.1) and (4.2). The transport calculation terminates when the maximum number of iterations or the maximum simulation time is reached. It is possible to restart from the final state of a previous T3D run, which is invaluable when a transport steady-state solution (characterised by power and particle losses that match sources and profiles that are constant in time) is not reached within a single run.

Figure 13.

Workflow of T3D-GS2. The main cycle is decomposed in four main blocks (in black): writing GS2 input files, running GS2, computing the quasi-linear heat flux and evolving the density and temperature profiles. The blocks in grey are executed only once per simulation.

4.2. Evolution of density and temperature profiles in a STEP case

Here, we present results from the first flux-driven transport simulations to account for hybrid-KBM turbulence in a STEP flat-top operating point.

4.2.2. Transport simulation results

Figure 14 compares the initial temperature, density and pressure profiles and their gradients taken from the JINTRAC STEP flat-top operating point with the corresponding profiles from the final state of the reference T3D simulation. Fluxes and gradients in T3D are computed at the mid-point between each radial grid point. The final electron temperature profile values are lower than in the initial condition, with the largest reduction occurring near the magnetic axis where $T_e$ decreases by approximately 10 %. The final normalised logarithmic electron temperature gradient, $a/L_{T_e}$, oscillates radially around its initial profile, and the largest difference is in the inner core, where the final electron temperature profile is flatter than the initial one.

Figure 14.

Comparison of initial (dotted lines) and final (solid lines) temperature (a), density (c) and pressure (e) profiles as well as their inverse gradient scale lengths (b,d,f). The initial profiles are taken from the JINTRAC analysis (Tholerus et al. Reference Tholerus, Casson, Marsden, Wilson, Brunetti, Fox, Freethy, Hender, Henderson and Hudoba2024). The round markers denote the position of the radial grid points in the T3D grid. The value of the density and pressure at the outermost radial grid point is imposed by the finite Dirichlet boundary condition. The shaded area represents the profile variation corresponding to a $\pm$40 % variation in $\gamma _E$ and a $\pm$20 % variation in $Q_0$ and $\alpha$. The solid line represents an interpolation through the T3D radial grid points. The particle source (star markers) as well as the electron (square markers) and ion (diamond markers) total power sources are shown in (c,e).

A significant reduction is observed in the ion temperature profile for $r/a<0.4$, where there is a 20 % reduction in $T_i$ near the magnetic axis. On the other hand, the initial and final ion temperature profiles are very similar in the outer core at $r/a>0.5$. The ITG is smaller than the initial temperature gradient over the entire radial domain, except near $r/a\simeq 0.8$, where it is slightly larger. We note that the initial ion temperature profile is significantly higher than the electron temperature despite this flat-top operating point being dominated by electron heating: this arises because the transport model used to design STEP-EC-HD assumed negligible turbulent ion heat transport. This feature is not recovered in the T3D calculation using our reduced model for hybrid-KBM turbulence, where the ion turbulent transport is non-negligible and the final $T_i$ and $T_e$ profiles are very similar, with $a/L_{T_i}$ reducing significantly compared with its initial value between $r/a\simeq 0.4$ and $r/a\simeq 0.65$.

The final T3D electron density is close to the initial profile for $r/a \leq 0.6$, but higher for $r/a>0.6$: i.e. there is a small increase in the edge density gradient for $r/a\simeq 0.8$ and a relaxation of the gradient around $r/a\simeq 0.65$. We note that the density increase occurs in proximity of the maximum density source due to pellet injection, shown in figure 14(c). The particle confinement time evaluated from the T3D steady-state solution is $\tau _p = \int n_e \,\mathrm {d} V / \int S_n\,\mathrm {d} n = 12.5$ s, which is close to that assumed in STEP-EC-HD (Tholerus et al. Reference Tholerus, Casson, Marsden, Wilson, Brunetti, Fox, Freethy, Hender, Henderson and Hudoba2024). We highlight that the T3D steady-state ion temperature and density gradients at $r/a\simeq 0.8$ are larger than in the initial condition from STEP-EC-HD. This may appear counter-intuitive, especially given that the reduced model heat and particle fluxes are much larger than the available sources at the STEP-EC-HD initial condition, as illustrated in figure 16. On the other hand, the dependence of the turbulent fluxes on the temperature and density gradients is non-monotonic in STEP-EC-HD because of the stabilising $\beta '$ effect, as pointed out in Giacomin et al. (Reference Giacomin, Kennedy, Casson, Ajay, Dickinson, Patel and Roach2024). The profile evolution in the initial phase leads to a steepening of temperature and density profile in the outer core region, which causes a local increase of the $\beta '$ value and a subsequent reduction of turbulent fluxes, thus sustaining larger density and temperature gradients in the outer core region. We note, however, that the variation of other parameters, such as $\beta$ and collisionality, the evolution of the fusion power and the contribution from the neoclassical fluxes make the dynamics regulating the profile evolution strongly nonlinear: this effect cannot be solely explained by the dependence of turbulent fluxes on gradients.

Figure 14 shows also a comparison between the initial and final pressure profile: the reduced final pressure for $r/a<0.6$ is consistent with lower $T_e$ and $T_i$ in this region of the plasma; for $r/a>0.6$ both the electron and ion final pressure profiles are larger than the initial profiles, mostly because of the increased density. If we do not separate transport and radiation losses, the energy confinement time using the final steady-state pressure profile gives $\tau _E = \frac {3}{2}\int p \,\mathrm {d} V / P_\mathrm {net} = 1.5$ s, with $P_\mathrm {net} = P_\mathrm {\alpha } + P_\mathrm {EC}$. If we exclude radiation losses from the net loss power using $P^{*}_\mathrm {net} = P_\mathrm {\alpha } + P_\mathrm {EC} - P_\mathrm {rad}$, this gives $\tau ^*_E = \frac {3}{2}\int p \,\mathrm {d} V / P^{*}_\mathrm {net} = 4.6$ s. These energy confinement times are similar to those assumed in the original STEP flat-top operating pointFootnote ¹² (Tholerus et al. Reference Tholerus, Casson, Marsden, Wilson, Brunetti, Fox, Freethy, Hender, Henderson and Hudoba2024). We highlight, however, that the radial domain evolved by these T3D simulations does not include the pedestal region. The outermost radial grid point boundary condition of the density and temperature is taken from near the pedestal top of the JINTRAC initial profiles (see figure 14); the pedestal model used in the JINTRAC modelling of this STEP flat-top operating point is described in Tholerus et al. (Reference Tholerus, Casson, Marsden, Wilson, Brunetti, Fox, Freethy, Hender, Henderson and Hudoba2024). Any variation of the pedestal top, and therefore of the outer boundary condition, is likely to impact core transport and the energy confinement time. Assessing the sensitivity on the outer boundary condition, i.e. on the density and temperature values at the pedestal top, is clearly a very important topic for future investigation.

The shaded area in figure 14 indicates the range of profiles variation obtained from T3D simulations with local values of $\gamma _E$ varied by $\pm$40 % from the nominal value, and also includes additional T3D simulations where $Q_0$ and $\alpha$ are varied by $\pm$20 % to assess the impacts of uncertainties in the reduced model parameters. The variation in the density and temperature profiles is comparatively small, which may seem surprising given the large sensitivity of fluxes to $\gamma _E$ in STEP-EC-HD. This relatively weak sensitivity of the kinetic profiles on the shearing rate can be explained (i) by T3D converging to a regime where the hybrid-KBM turbulence is closer to marginal, especially at low $k_y$, in the transport steady state, owing reductions of the density gradient in the core and increases in $\beta '$ at the edge; and (ii) by large transport stiffness, where small variations in the temperature and density gradients lead to large changes in the turbulent flux. High transport stiffness improves the robustness of the T3D steady state.

Figure 15 shows the time evolution of the density and temperature gradients as well as of the power and particle losses over the duration of the T3D reference simulation. The time traces can be divided in three main phases. The first phase is very short ($0< t\lesssim 0.1$ s, 50 solver iterations) and is characterised by a fast evolution of the profiles that leads to a substantial reduction of power and particle losses. This fast initial phase is caused by the strong initial particle and power imbalance, as shown in figures 16(a) and 16(c). The second phase ($0.1\lesssim t\lesssim 9$ s, 100 solver iterations) is characterised by a slow evolution of gradients and fluxes towards a steady-state solution. In this phase, the particle fluxes at $r/a\simeq 0.16$ and $r/a\simeq 0.33$ drop to negligible values, which is consistent with an approach to steady state with no particle source in the core. There is a transitory ripple in the time evolution of the power and particle losses at $r/a\simeq 0.33$ and $r/a\simeq 0.49$, while the density and temperature gradients vary smoothly; this is caused by low $k_y$ hybrid-KBMs being stabilised or destabilised while T3D profiles are slowly converging to a steady-state solution, and it is amplified by stiff transport. Apart from the ripple and slow increase of particle flux at $r/a\simeq 0.5$, power and particle losses do not vary significantly in the second phase. The final phase (40 solver iterations) reaches the transport steady state, where the kinetic gradients, power and particle losses are constant in time. Figure 16 shows excellent agreement between transport losses and sources in the T3D transport steady state.

Figure 15.

Time evolution of the density and temperature gradients as well as of the power and particle losses in the reference T3D simulation. The initial profiles are taken from the JETTO STEP-EC-HD flat-top operating point. Different colours represent different radial surfaces (fluxes and gradients in T3D are computed at the mid-point between each radial grid point).

Figure 16.

Power (a,b) and particle (c,d) balance in the initial and final states of the reference T3D simulation. Heat and particle sources are denoted with open markers, and heat and particle losses with solid markers. The initial state is taken from JINTRAC-JETTO STEP flat top.

The heating power density radial profiles from auxiliary and $\alpha$-particle heating, collisional exchange and radiation are shown in figure 17(a) for the final T3D steady state, while the evolution of the fusion power is shown in figure 17(b). Consistently with the pressure profile evolution, the fusion power rapidly decreases in the initial phase, which is followed by a slower evolution. The final fusion power is approximately 10 % smaller than its initial value.

Figure 17.

(a) Power source radial profile in the final T3D steady state of the auxiliary heating ($S_\mathrm {aux}$), $\alpha$-particle heating ($S_{\alpha, e}$ and $S_{\alpha, i}$), collisional exchange ($S_{\mathrm {coll}, e}=-S_{\mathrm {coll}, i}$) and radiation ($S_\mathrm {rad}$). (b) Fusion power evolution in the reference T3D simulation. The fusion power is normalised to its initial value of 1.7 GW.

5.1. Linear analysis

Linear simulations are performed on these three surfaces with GS2, using numerical resolutions similar to those used in § 2.1, except setting $n_\text {period} = 32$ at $r/a=0.5$. The dominant mode growth rates and mode frequencies are shown as functions of $k_y$ in figure 18. The mode frequency at $r/a=0.5$ shows a discontinuity between $k_y\rho _s=0.15$ and $k_y\rho _s=0.2$, which corresponds to a transition of the dominant mode from hybrid-KBMs to MTMs. The $\delta \phi$ parallel mode structures at $k_y\rho _s \simeq 0.07$ and $k_y\rho _s \simeq 0.83$, corresponding to the fastest growing hybrid-KBM and MTM respectively, are shown in figure 19. The hybrid-KBM eigenfunction is narrow in $\theta$ and has even parity, while the MTM eigenfunction has odd parity and is much more extended in $\theta$. A sign change on the mode frequency is also observed in figure 18 at $r/a=0.8$ around $k_y\rho _s \simeq 0.2$. In this case, however, the parallel structures of modes with positive and negative frequencies are similar. This is consistent with the complex nature of the hybrid-KBM, involving both trapped electron and ion dynamics, discussed in Kennedy et al. (Reference Kennedy, Giacomin, Casson, Dickinson, Hornsby, Patel and Roach2023).

Figure 18.

Growth rate (a) and mode frequency (b) as functions of $k_y$ at $r/a=0.5$ (blue dots), $r/a=0.65$ (orange dots) and $r/a=0.8$ (green dots) from linear simulations of the new STEP density and temperature profiles.

Figure 19.

Parallel mode structure of $\delta \phi$ at $k_y\rho _s\simeq 0.07$ (a) and $k_y\rho _s\simeq 0.83$ (b) corresponding to the maximum growth rate of the hybrid-KBM and MTM instability, respectively, at $r/a=0.5$ of the new STEP profiles obtained from the T3D reference simulation.

Figure 20 shows a comparison of the growth rate values between the initial and final STEP plasma profiles at $r/a=0.5$, $r/a=0.65$ and $r/a=0.8$. Beyond the difference on the dominant instability at $r/a=0.5$ (the hybrid-KBM and MTM instabilities dominate in the original STEP case and in the new T3D steady state, respectively), a strong stabilisation at low $k_y$ and a substantial growth rate maximum shift at higher $k_y$ are observed in the new T3D steady state. At both $r/a=0.65$ and $r/a=0.8$ cases, where the dominant instability remains the hybrid-KBM, there is a clear reduction of the growth rates across most of the spectrum, including at low $k_y$ which typically dominates the turbulent fluxes. This linear comparison underlies the substantial reduction of turbulent flux predictions by the reduced transport model between the initial condition and the transport steady state. This will now be compared with nonlinear gyrokinetic simulations.

Figure 20.

Growth rate comparison between the initial (STEP-EC-HD) and final (T3D-GS2) profiles at $r/a=0.5$ (a), $r/a=0.65$ (b) and $r/a=0.8$ (c).

5.2. Nonlinear analysis

Here, we briefly show results from nonlinear gyrokinetic simulations performed using GENE at $r/a = [0.5, 0.65, 0.8]$ in the T3D steady-state solution: these simulations used the numerical resolutions listed in table 3. These simulations include equilibrium flow shear with $\gamma _E = [0.02,0.03,0.05 ]\,c_s/a$ at $r/a = [0.5,0.65,0.8]$, respectively. Figure 21 plots the total heat and particle transport losses across these three surfaces in the transport steady state, comparing values from the quasi-linear model in T3D with neoclassical and turbulent fluxes computed using NEO and GENE.

Figure 21.

Turbulent (red dots), neoclassical (blue dots) and total (green dots) power (a) and particle (b) losses at $r/a=0.5$, $r/a=0.65$ and $r/a=0.8$ of the new STEP profiles, compared with the expected heating power and particle fuelling from the original STEP case (black dots). The green shaded area corresponds to a $\pm$20 % of the total heat (a) and particle (b) transport. The turbulent fluxes are evaluated from GENE nonlinear gyrokinetic simulations and the neoclassical fluxes are evaluated from NEO simulations, both considering the T3D steady-state profiles.

Table 3.

Numerical resolutions used in nonlinear GENE simulations on three surfaces from the T3D transport steady state, $r/a=[0.5, 0.65,0.8]$.

At $r/a=0.5$ and $r/a=0.65$, figure 21 shows relatively good agreement between the nonlinear gyrokinetic heat flux and the reduced transport model in T3D. The neoclassical heat flux is negligible at these radial locations. There is also good agreement between the gyrokinetic and reduced model particle losses here, and we note that both neoclassical and turbulent particle fluxes are negligible at $r/a=0.5$.

Turbulent transport in the GENE nonlinear simulation at $r/a=0.5$ is mainly driven by hybrid-KBMs, in spite of linearly dominant MTMs (see figure 18): the nonlinearly dominant mode fields $\delta \phi (\theta )$ and $\delta A_\parallel (\theta )$ have even and odd parity, respectively. Good agreement in figure 21 at $r/a=0.5$ suggests that transport from these hybrid-KBMs is captured faithfully by the quasi-linear reduced transport model.

At $r/a=0.8$, the quasi-linear model overestimates the heat and particle fluxes from the nonlinear simulation, despite turbulent transport being solely driven by hybrid-KBMs. Turbulent fluxes should therefore be described well by the quasi-linear model. The disagreement might be due to some inaccuracy affecting the quasi-linear model, as visible in figure 10 where a few points exhibit order of magnitude differences between the quasi-linear and the nonlinear heat flux predictions. In addition, the local equilibrium at $r/a=0.8$ might be near a marginal state where a small increase in drive can boost the turbulent fluxes significantly: this marginal condition, where simulations are more difficult and nonlinear effects intrinsically important, might not be described well by the quasi-linear model, although further investigations, which are outside the aim of this work, would be required before drawing firm conclusions. We note that transport stiffness should allow particle and power balance to be restored with a modest increase in the drive and minimal impact on the overall density and temperature profiles. In fact, transport stiffness makes flux-driven simulations less sensitive to input parameters and improves the robustness of the steady-state profile predictions.

In summary, figure 21 shows that gyrokinetic transport fluxes for the T3D transport steady state in STEP are broadly consistent with the quasi-linear reduced transport model of (3.7) and (3.8). This gives confidence in the validity of the reduced model within the parameter space explored in this T3D transport simulation. It is important to stress, however, that this is a first iteration to develop a suitable reduced model for STEP: further guidance from higher-fidelity nonlinear gyrokinetic simulations will be required to overcome important limitations of the reduced transport model presented here, including the absence of fast $\alpha$ effects. Nevertheless the reduced model and T3D-GS2 are powerful tools that could be further developed and used more extensively, together with existing integrated modelling codes, to assist in the design of new STEP plasmas and flat-top operating points.