Electric field and turbulence in global Braginskii simulations across the ASDEX Upgrade edge and scrape-off layer

In the confined region, to lowest order, we can assume $\left\langle \mathbf{v}_E \cdot \nabla n \right\rangle = 0$ since for the background we have only a radial density gradient and no radial E × B advection. For now, let us also neglect any flows, i.e. $\left\langle \nabla\cdot n(\mathbf{u}_\mathrm{pol} + u_\parallel\mathbf{b}) \right\rangle = 0$. In absence of stationary sources S_n , the balance for the stationary fields becomes

$$\begin{equation} n\nabla \cdot ( \mathbf{v}_E ) + \nabla \cdot ( n\mathbf{v}_{*}^\mathrm{i} ) = - nC(\varphi) - C(p_\mathrm{i}) / e = 0. \end{equation} \tag{ 9 }$$

As there are no background poloidal gradients along closed flux surfaces, $\partial_\theta f \approx 0$, the balance is fulfilled by the static equilibrium radial electric field

$$\begin{equation} E_r = \frac{\partial_r p_\mathrm{i}}{en}. \end{equation} \tag{ 10 }$$

We expect the dominant contribution to the radial electric field on closed flux surfaces to be this balance between E × B and diamagnetic compression of ion density along the curved magnetic field [1, 2]. Corrections, especially in L-mode, can be significant though – see section 4.1. It is worth pointing out that under this equilibrium electric field, the background plasma is at rest (i.e. static), but fluctuations are E × B advected. Deviations from equation (9), on the other hand, also impact the background.

Besides the lowest order equilibrium electric field (10), additional contributions arise from polarisation and parallel particle fluxes contained in equation (8). Obtaining closed form analytical expressions for their stationary contributions is not possible. In this section we discuss the form and role of the ion polarisation velocity and viscosity, while numerical simulation results will be presented in section 4.

As noted at equation (2), the polarisation velocity is the correction to the perpendicular ion velocity which contains inertia and viscosity. However, the latter depend themselves on the perpendicular ion velocity in a non-trivial way. Therefore, in reduced fluid models, the polarisation velocity is approximated to first order by inserting the zero order approximation $\mathbf{v}_{\perp0}^\mathrm{i} = \mathbf{v}_E+\mathbf{v}_{*}^\mathrm{i}$ into inertia and viscous stress. Since the divergence of the zeroth order velocities is of the same order as the divergence of the polarisation velocity, $\nabla\cdot\left(n\mathbf{u}_\mathrm{pol}\right) \sim \nabla\cdot\left(n\mathbf{v}_{\perp0}^\mathrm{i}\right)$, it has to be retained for divergence terms [36, 37]. Explicitly, the expression reads

$$\begin{align} \nabla \cdot \left( n \mathbf{u}_\mathrm{pol} \right) & = \nabla \cdot \left[ - \frac{n M_\mathrm{i}}{eB^2} \left( \frac{\partial}{\partial t} + \mathbf{v}_E \cdot \nabla + u_\parallel \nabla_\parallel \right)\right.\nonumber\\ & \quad \times\left. \left( \nabla_\perp \varphi + \frac{\nabla_\perp p_\mathrm{i}}{en} \right) \right] - \frac{1}{6e} C(G), \end{align} \tag{ 11 }$$

where the first term on the right hand side represents inertia, and the second viscosity. We can also write

$$\begin{equation} \nabla \cdot \left( n \mathbf{u}_\mathrm{pol} \right) = \nabla \cdot \left( n \mathbf{u}_\mathrm{pol}^\mathrm{in} \right) - \frac{1}{6e} C(G). \end{equation} \tag{ 12 }$$

The viscous stress function [38] is

$$\begin{equation} G = - \eta^\mathrm{i}_0 \left[ \frac{2}{B^{3/2}} \nabla \cdot \left( u_\parallel B^{3/2} \mathbf{b} \right) - \frac{1}{2} \left( C(\varphi) + \frac{1}{en}C(p_\mathrm{i}) \right) \right], \end{equation} \tag{ 13 }$$

with the dominant Braginskii viscosity coefficient $\eta^\mathrm{i}_0 = 0.96 p_\mathrm{i} \tau_\mathrm{i}$ and ion collision time $\tau_\mathrm{i}$.

Due to the time derivative in the inertial part, the direct implementation of the ion continuity equation is cumbersome. Instead, in drift-reduced Braginskii codes like GRILLIX, the electron continuity equation is used – where the polarisation velocity can be neglected due to the small electron mass – together with the quasi-neutrality, or vorticity, equation $e\partial_t(n_\mathrm{e}-n_\mathrm{i}) = \nabla \cdot \mathbf{j} = 0$. The latter reads

$$\begin{equation} e\nabla \cdot \left( n \mathbf{u}_\mathrm{pol} \right) = C( p_\mathrm{e} + p_\mathrm{i} ) - \nabla \cdot ( j_\parallel \mathbf{b} ), \end{equation} \tag{ 14 }$$

with the parallel current $j_\parallel = en ( u_\parallel - v_\parallel )$ and parallel electron velocity $v_\parallel$. If the right hand side of (14) balances, i.e. the diamagnetic charge separation is balanced by a parallel current, we obtain the Pfirsch–Schlüter current [39, chapter 8.4]. However, we will see in section 4 that the divergence of the polarisation flux is not necessarily small.

Equation (14) is typically called vorticity equation because under some (too) strong approximations, regarding the size of density fluctuations and homogeneity of the magnetic field, one can rewrite

$$\begin{equation} \nabla \cdot \left[ \frac{n}{B^2} \frac{\mathrm{d}}{\mathrm{d} t} \left( \nabla_\perp \varphi + \frac{\nabla_\perp p_\mathrm{i}}{en} \right) \right] \sim \frac{n_0}{B} \frac{\mathrm{d}}{\mathrm{d} t} \mathbf{b} \cdot \nabla \times \left( \mathbf{v}_E + \mathbf{v}_{*}^\mathrm{i} \right). \end{equation} \tag{ 15 }$$

The quantity on the right hand side is the 2D fluid vorticity $\Omega_\mathrm{F} = \mathbf{b} \cdot \nabla \times \left( \mathbf{v}_E + \mathbf{v}_{*}^\mathrm{i} \right)$, not to be confused with the generalised vorticity defined in appendix A. Note that while the diamagnetic velocity is 'nearly divergence-free', it is not rotation-free.

The inertial part of the vorticity equation is central to fluid plasma turbulence models, particularly the non-linearity $\mathbf{v}_E \cdot \nabla \mathbf{E}_\perp$, as it is ultimately responsible for the turbulent inverse energy cascade [40]. This term also contains the Reynolds stress, responsible for zonal flow drive [5]. We note that a rigorous decomposition into mean field and fluctuating flows in a global model is quite complex, e.g. one requires the Favre instead of the Reynolds average [41]. Further, a second Reynolds stress term contained in $u_\parallel \nabla_\parallel \nabla_\perp \varphi$ couples 2D turbulence to parallel dynamics [42]. This decomposition will not be further detailed here. In the steady-state, we must have $\left\langle \partial_t \left( \nabla_\perp \varphi + \frac{\nabla_\perp p_\mathrm{i}}{en} \right) \right\rangle = 0$. Therefore, we will use $\left\langle \nabla \cdot \left( n \mathbf{u}_\mathrm{pol}^\mathrm{in} \right) \right\rangle$ as a measure of the stationary Reynolds stress.

In this section we summarize some key insights from neoclassical theory [39] for the confined plasma. The main difference to our model is that inertia, and hence the polarisation velocity $\mathbf{u}_\mathrm{pol}$ from the previous chapter, is neglected. In this case it is then a good assumption that pressure and electrostatic potential are constant along closed flux surfaces ψ. With $\partial_r f = RB_p \partial_\psi f$, one obtains the relation [39, chapter 8.5]

$$\begin{equation} E_r = \frac{1}{en}\frac{\partial p_\mathrm{i}}{\partial r} + u_\parallel B_\theta \frac{B}{B_\phi} - v_\theta \frac{B^2}{B_\phi}. \end{equation} \tag{ 16 }$$

$v_\theta$ is the flux surface averaged poloidal rotation. An important result is then that according to the Braginskii fluid closure, poloidal rotation is damped to zero due to the ion viscosity [39, chapter 12.3]. The reason is that in a flux surface average over the parallel momentum equation (A3), one finds the condition

$$\begin{equation} \left\langle \mathbf{B} \cdot \nabla \cdot \Pi_\mathrm{i} \right\rangle_S = 3\eta^\mathrm{i}_0 \left\langle \left( \nabla_\parallel B \right)^2 \right\rangle_S v_\theta = 0. \end{equation} \tag{ 17 }$$

$\Pi_\mathrm{i}$ is thereby the Braginskii ion viscous stress tensor. We will write this with the viscous stress function G defined above in equation (13),

$$\begin{equation} \mathbf{B} \cdot \nabla \cdot \Pi_\mathrm{i} = \frac{2}{3}B^{5/2}\nabla_\parallel \frac{G}{B^{3/2}}. \end{equation} \tag{ 18 }$$

The situation is more complicated, however, when turbulent transport is considered (i.e. $\mathbf{u}_\mathrm{pol}$) – see also [39, chapter 13.1] – as will be detailed in section 4.1 on the basis of numerical simulations. Note, as the viscous stress function G appears in Braginskii's ion temperature equation (A6), this flow damping acts as a heating mechanism for ions known as 'magnetic pumping'.

In the SOL, due to the absence of closed flux surfaces, parallel dynamics and boundary conditions take a more important role [9, 10]. The electrostatic Ohm's law – used here for simplicity, although our full model is electromagnetic, see appendix A – reads

$$\begin{equation} \eta_\parallel j_\parallel = - \nabla_\parallel \varphi + \frac{\nabla_\parallel p_\mathrm{e}}{en} + \frac{0.71}{e} \nabla_\parallel T_\mathrm{e}, \end{equation} \tag{ 19 }$$

with resistivity $\eta_\parallel$. We have insulating sheath boundary conditions $j_\parallel = 0$ at the divertor. And in present simulations, without any neutral particles or impurities, there are no significant sources or sinks for any terms in the above equation. This results in flat parallel gradients (in an average over turbulent fluctuations). Therefore, $T_\mathrm{e}$ is roughly equal at the outboard mid-plane and at the divertor, and resistivity is negligibly small. At the divertor, the insulating sheath boundary condition

$$\begin{equation} \varphi = \Lambda T_\mathrm{e} \end{equation} \tag{ 20 }$$

for the potential applies, with $\Lambda = -\frac{1}{2} \ln \left[\! \left(\, 2\pi\frac{m_\mathrm{e}}{M_\mathrm{i}} \!\right) \left(\! 1 + \frac{T_\mathrm{i}}{T_\mathrm{e}} \!\right)\!\right] \approx 2.69$. Hence, due to $\nabla_\parallel \varphi \approx 0$, we expect $\varphi = \Lambda T_\mathrm{e}$ to also hold at the outboard mid-plane. Having $E_r \sim \partial_r p_\mathrm{i}/en$ in the confined region (10) and $E_r \sim -\Lambda \partial_r T_\mathrm{e}$ in the SOL, we expect a jump of the electric field across the separatrix, and contra-rotating poloidal flows inside and outside of it [11].

It should be noted that in reality, contrary to our yet simplified model, the presence of neutral gas and impurities in the SOL leads to more complex parallel gradients, even in attached conditions. Furthermore, the sheath can conduct significant currents. Therefore, the SOL dynamics itself can be much more involved, particularly in detached plasmas [43].

This section details the setup of the simulations – namely geometry, parameters, initial state, resolution and computational cost – guided by the AUG discharge #36190 at time t = 2–4 s. In this discharge the toroidal magnetic field was $B_\mathrm{tor} = -2.5$ T, i.e. in favourable configuration with $\mathbf{B}\times\nabla B \sim - \hat{\mathbf{e}}_Z$ towards the X-point. The plasma current was I_p  = 0.8 MA ⁵ resulting in q₉₅ = 4.4, and the average triangularity was δ = 0.21. Further important inputs are the measured $T_\mathrm{e} \approx 350\,\,{\textrm{eV}}$ and $n_\mathrm{e} \approx 2 \times 10^{19}$ m⁻³ at $\rho_\mathrm{pol} = 0.9$, and $T_\mathrm{e} \approx 50$–80 eV and $n_\mathrm{e} \approx 1 \times 10^{19}$ m⁻³ at the separatrix. The total heating power in the discharge was roughly 800 kW: 550 kW neutral beam injection, 500 kW ohmic heating and subtracting 250 kW radiation losses.

The geometry is illustrated in figure 1 by the poloidal flux function Ψ. The separatrix and machine walls are highlighted. The simulation domain extends beyond the divertor legs due to the implementation of parallel boundary conditions via penalisation [33]. We define the normalized poloidal flux radius as $\rho_\mathrm{pol} = \sqrt{\frac{\Psi - \Psi_0}{\Psi_X - \Psi_0}}$, with Ψ₀ and Ψ_X the poloidal magnetic flux at magnetic axis and separatrix, respectively. Simulations are performed for $\rho_\mathrm{pol} \gt 0.9$, since our fluid model is not valid in the plasma core, and the reduced domain saves computational cost. The outer most flux surface is chosen at the location where the HFS main chamber wall acts as a plasma limiter, at $\rho_\mathrm{pol} = 1.05$.

Zoom In Zoom Out Reset image size

Our choice of reference values is guided by parameters and measurements in AUG discharge #36190: major radius R₀ = 1.65 m, magnetic field on axis B₀ = 2.5 T, density $n_0 = 10^{19}$ m⁻³, electron and ion temperatures $T_\mathrm{e0} = T_\mathrm{i0} = 100$ eV, $Z_\mathrm{eff} = 1.3$ and deuterium ion mass $M_\mathrm{i} = 2m_\mathrm{p}$. We note that in a global model, reference density and temperature are only required for the normalization of the equations, to write them in dimensionless form as in appendix A, and have no physical relevance. They are chosen only on the order of measured separatrix values, since their exact choice does not matter. $Z_\mathrm{eff}$ is considered only in the calculation of the electron collision frequency. The resulting dimensionless collisionless parameters of the system – see appendix A for their definition – are δ = 2854.2, $\beta_0 = 3.227 \times 10^{-5}$, µ = 2.723 × 10⁻⁴ and ζ = 1. The collisional parameters are $\nu_\mathrm{e0} = 12.30$, $\eta_{\parallel 0} = 0.0017$, $\chi_{\parallel \mathrm{e}0} = 940$, $\chi_{\parallel \mathrm{i}0} = 35.35$ and $\eta_\mathrm{i0} = 8.70$. The Braginskii heat conductivity is known to be inappropriate at low collisionality [22, 44, 45]. For instance, it becomes arbitrarily large as ∼T^5/2 due to the lack of kinetic flux limiting effects, such as Landau damping. In this work, we use a simple flux limiter which reduces the stiffness and therefore computational expense of the parallel heat conduction: $\chi_{\parallel \mathrm{e,i} 0} T_\mathrm{e,i}^{\,5/2} \leqslant 940$. For the chosen reference values, we obtain as reference drift scale and ion Larmor radius

$$\begin{equation} \rho_{s0} = \frac{\sqrt{M_\mathrm{i} T_\mathrm{e0}}}{eB_0} = \rho_\mathrm{i0} = \frac{\sqrt{M_\mathrm{i} T_\mathrm{i0}}}{eB_0} = 0.578\,\textrm{mm}. \end{equation} \tag{ 21 }$$

Note that the local Larmor radius varies with magnetic field and temperature. E.g. at reference temperature, it is larger by 25% at the outboard mid-plane due to the reduced local magnetic field of B ≈ 2 T.

In diverted geometry and in absence of neutral gas ionization, we find that density can drop arbitrarily low in the far SOL and private flux region. This contradicts the experimental observation that density can even rise in the far SOL [46], i.e. the simulations are missing the necessary mechanism, e.g. neutral gas ionization. It also makes the solution of the equations unnecessarily expensive, since in particular the parabolic part becomes very stiff at low collisionality. As a simple solution we apply an adaptive source to keep the density above 3 × 10¹⁷ m⁻³, and temperatures above 3 eV. This source does not hinder cross-field inflow but limits parallel outflow in the far SOL, starting to act above $\rho_\mathrm{pol} = 1.02$ in the reference simulation, but e.g. only above $\rho_\mathrm{pol} = 1.04$ in the simulation in section 5.2. For simplicity, respectively, due to the unphysical core boundary at $\rho_\mathrm{pol} = 0.9$, simulations are adaptively flux driven: the heat and particle sources at the core boundary are adapted in time such as to hold density and temperature fixed there – while in the rest of the domain they evolve freely.

In a global code, background profiles are not chosen but rather evolve freely in accordance with turbulent fluxes. Nonetheless, an initial state must be chosen for the simulation. The most trivial choice would be a flat, or even zero, background profile – which then evolves due to sources at the core boundary. However, in our experience the saturated state is reached faster the closer the initial state is to the final. We choose as initial state for density and temperatures the simple profile as function of $\rho_\mathrm{pol}$

$$\begin{align} f(\rho) = \begin{cases} f_\mathrm{ped} & \textrm{for: } \rho \lt \rho_\mathrm{ped}, \\ -A\sin(a\cdot \rho + b) + B & \textrm{for: } \rho_\mathrm{ped} \leqslant \rho \leqslant \rho_\mathrm{sep}, \\ f_\mathrm{sep} & \textrm{for: } \rho_\mathrm{sep} \lt \rho, \end{cases} \end{align} \tag{ 22 }$$

with

$$\begin{align*} A = & \,\,\frac{f_\mathrm{ped} - f_\mathrm{sep}}{2}, \qquad B = \,\, f_\mathrm{ped} - A, \\ a = \,\, & \frac{\pi}{\rho_\mathrm{sep} - \rho_\mathrm{ped}},\quad b = -\frac{\pi}{2}\frac{\rho_\mathrm{sep}+\rho_\mathrm{ped}}{\rho_\mathrm{sep}-\rho_\mathrm{ped}}.\\[-6pt] \end{align*}$$

In present simulations, we chose $\rho_\mathrm{ped} = 0.92$ and $\rho_\mathrm{sep} = 0.999$. The value $f_\mathrm{ped}$ is held constant between ρ = 0.9 and ρ = 0.92 by an adaptive source. In the SOL, the initial profile is flat and equal to $f_\mathrm{sep}$.

The core boundary, respectively, pedestal top values are chosen in accordance to the experiment as $n_\mathrm{ped} = 2 \times 10^{19}$ m⁻³ and $T_\mathrm{ped}^\mathrm{\,e,i} = 350$ eV. The separatrix and SOL initial values are less relevant because they adapt in the course of the simulation. An important restriction comes from the choice of flat SOL profiles: simulations saturate faster with $f_\mathrm{sep}$ as low as possible, since it also initialises lower far SOL values. On the other hand, ideal ballooning stability of the initial profiles prohibits large radial gradients in the confined region, which is correlated with the Greenwald density limit [47]. Our choice for the reference case is $n_\mathrm{sep} = 4 \times 10^{18} \,{ }\textrm{m}^{-3}$, $T_\mathrm{sep}^\mathrm{\,e} = 40$ eV and $T_\mathrm{sep}^\mathrm{\,i} = 50$ eV. In the confined region, electrostatic potential and current are chosen such that the initial pressure profile is stable: $\partial_r \varphi = - \partial_r p_\mathrm{i} / (en)$ according to (10), resulting in $\mathbf{u}_\mathrm{pol} = G = \Omega = 0$ (no vorticity and polarisation velocity), and $\nabla \cdot ( j_\parallel \mathbf{b} ) = eC( p_\mathrm{e} + p_\mathrm{i} )$ according to (14). The parallel velocity is chosen as zero. In the SOL, the initial state is chosen such that it fulfils the Bohm boundary conditions (A10). Particularly, the parallel velocity is $u_\parallel = \pm \sqrt{(T_\mathrm{e} + T_\mathrm{i})/M_\mathrm{i}}$ at each divertor plate, with a linear dependence on distance between the plates. The initial outboard mid-plane profiles are illustrated in figure 2. To trigger an instability and consequent turbulence, random noise of magnitude 10⁻⁵ is added to the density and temperature profiles.

Zoom In Zoom Out Reset image size

For the reference simulation, we choose as poloidal resolution $h_f = 2.5 \rho_{s0} = 1.45$ mm. Toroidally, we resolve 16 planes. The time step is chosen slightly below the stability limit as $\Delta t = 5 \times 10^{-5} R_0 / c_{s0} = 1.2$ ns. Due to the limited resolution, the turbulent spectrum must be cut at the poloidal grid scale. To this end, third order hyperviscosity is applied as detailed in appendix A, with $\nu_\perp = 3350$ for all fields. Being able to resolve even smaller scales is highly desirable, particularly in a global code. Eventually, our choice of resolution is a compromise, restricted by the allowable computational expense. Nevertheless, we also performed convergence tests with $h_f = 1.67\rho_{s0}$, and a separate test with 24 planes at 2.5ρ_s0. For the 1.67ρ_s0 simulation, hyperviscosity was reduced to 300. For the 24 planes simulation, the timestep had to be reduced to 0.7 ns.

It is worth mentioning the computational cost of these simulations, and how it scales. The reference case with 2.5ρ_s0 poloidal resolution and 16 toroidal planes (∼7 million points) required for simulating 1 ms of the discharge 9 days on 384 processor cores (8 nodes) of the Marconi-A3 SKL partition. As we find saturation only after $t \gtrsim 5$ ms, at least roughly 0.4 MCPUh and 2 months of patience were required per simulation. Increasing poloidal resolution to 1.67ρ_s0 raised both the number of grid points and the cost and duration of the simulation by a factor of 2. Increasing toroidal resolution raises the cost quadratically with the number of planes as the maximum allowed timestep is inversly proportional to toroidal resolution, i.e. 24 instead of 16 poloidal planes roughly doubles the cost.

We find no large differences in heat transport with increasing resolution, suggesting that it is mostly captured by the reference case. In detail, however, ion heat transport slightly increases with toroidal resolution, and zonal flows increase with poloidal resolution, as detailed in sections 3.2 and 4, respectively.

In current simulations, fixed density and temperature are prescribed at an artificial core boundary at $\rho_\mathrm{pol} = 0.9$ – held constant by an adaptive heat and particle source. In figure 3 we show how the input power $\left\langle P_\mathrm{e,i} \right\rangle_V = 3/2\left\langle T_\mathrm{e,i} S_n + n S_{T_\mathrm{e,i}} \right\rangle_V$ for the reference simulation varies in time – with particle and temperature sources S_n and $S_{T_\mathrm{e,i}}$ averaged over the whole volume. In the same plot, we show the globally averaged electrostatic potential oscillations (in units of 25 V) – this is due to the GAM, which is always present in global edge simulations [27]. It is important that flux-surface-symmetric potential oscillations like the GAM are permitted up to the core boundary through the implementation of the zonal Neumann boundary condition (see appendix A). But the GAM also leads to an up-down asymmetric pressure oscillation [48], while the adaptive source acts to hold pressure constant – and therefore follows the GAM with the same frequency, damping the pressure oscillation at the core boundary. While the oscillation of the source is a numerical artifact, on average $\left\langle P_\mathrm{e,i} \right\rangle_V$ saturates at $t \gtrsim 5$ ms on a reasonable level: 192 ± 29 kW for electrons and 224 ± 36 kW for ions, with one standard deviation quantifying the oscillations. This is roughly half the input power in AUG discharge #36190.

Zoom In Zoom Out Reset image size

In figure 4, we show zonally averaged separatrix values for density, electron and ion temperatures and electrostatic potential. As indicated by the input power, in the confined region turbulence, and largely also profiles, saturate within 3 ms of simulation time. Further outside, however, time scales are longer due to lower temperatures – particularly in the SOL, where the transit time is of the order 1 ms (as can be estimated from ∼ 50 m connection length and 60 km s⁻¹ flow velocity from figure 5). Turbulence at this time is locally saturated, but globally profiles vary predominantly due to the outflow in the SOL, as we started from flat profiles there. Therefore, for the confined region, the separatrix values are a good criterion for saturation. In the following, we will use the data from $t\gt4$ ms for our statistical analysis, although temperatures are still very slowly evolving. Additionally, the globally averaged parallel velocity is shown which is a measure of toroidal rotation. It saturates at 13 km s⁻¹ in the direction of the plasma current – a reasonable value for AUG [49].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In unfavourable configuration, i.e. with $\mathbf{B}\times\nabla B \sim + \hat{\mathbf{e}}_Z$ away from the X-point, the input power is slightly lower: 187 ± 29 and 204 ± 33 kW for electrons and ions, respectively. The simulation with 24 toroidal planes (higher toroidal resolution) has 176 ± 26 kW for electrons and 221 ± 32 kW for ions, respectively. For the 1.67ρ_s0 simulation (higher poloidal resolution), we get 187 ± 13 kW for electrons and 221 ± 17 kW for ions – i.e. no significant difference to the reference case with 2.5ρ_s0. This is remarkable as zonal flows (a radial modulation of the electric field, pressure and parallel velocity) are much more pronounced in the 1.67ρ_s0 simulation, as discussed in section 4, suggesting that overall thermal transport seems to be barely affected by zonal flows. The reduced oscillation amplitude also indicates a less pronounced GAM. There are no other qualitative differences between the reference simulation and the simulations with higher poloidal or toroidal resolution, as profiles and fluctuation levels remain very similar. Therefore, we conclude that thermal transport is largely resolved in the reference simulation. Increasing resolution further is desirable, but requires a significant speed-up of the code.

Figure 5(a) shows a snapshot of the plasma density in a poloidal plane, in quasi-steady state at $t = 308 R_0 / c_{s0} \approx 7$ ms for the reference simulation. Figure 5(b) displays profiles at the outboard mid-plane (Z = 0), averaged toroidally and in time over 1 ms, for density, temperature and parallel velocity. While the density profile is rather flat, we see steep pedestals in temperatures around the separatrix – the reason is high heat conductivity, as discussed in section 5.2. Experimentally, pedestals are observed in L-mode at low density, particularly for $T_\mathrm{e}$ [50].

We see some penetration of parallel velocity into the confined region, much less than without ion viscosity in figure 11 though. 'w_GTi  = 1' stands for active viscous ion heating (in equation (A6)) in the confined region – ions are significantly hotter in that case. It was not possible to run stable simulations with this effect active also in the SOL, hence why the reference simulation is without viscous ion heating – this is further discussed in section 5.1. If active also in the SOL, the heating is even stronger there, leading to non-monotonous $T_\mathrm{i}$ profiles peaking in the SOL and unstable simulations.

Figure 5(c) shows fluctuation levels at the outboard mid-plane. The potential fluctuations are normalised to background (mean) electron temperature, and are the largest in the system. As we have $\sigma_\varphi / \bar{T}_\mathrm{e} \gtrsim 2 \sigma_n / \bar{n}$ in the confined region, with $T_\mathrm{e}$ and $T_\mathrm{i}$ fluctuations an order of magnitude smaller, this indicates that turbulence in the confined region is ballooning driven [51, 52]. Around the separatrix, there is a local minimum in $\sigma_n / \bar{n}$ and $\sigma_\varphi / \bar{\varphi}$, indicating some level of turbulence suppression due to the sheared E × B flow. In the SOL, fluctuation amplitudes peak around ${\rho _{{\text{pol}}}} \approx 1.01{\text{-}}1.02$, at the bottom of the steep gradients. We have here $\sigma_\varphi / \bar{T}_\mathrm{e} \gg \sigma_n / \bar{n}$, roughly an order of magnitude, but also rather high $\sigma_{T_\mathrm{i}} / \bar{T}_\mathrm{i}$ and a steep ion temperature gradient. This suggests a combination of the Kelvin–Helmholtz instability [53, table I] and ITG drive.

The induction of a perturbed magnetic field, $\beta_0\partial_tA_\parallel$, is crucial for the dynamics of the parallel electric field and current, and therefore electrostatic turbulent transport. The perturbation $\tilde{\mathbf{B}} = \nabla\times A_\parallel\mathbf{b}$ does not necessarily lead to significant additional electromagnetic transport, though: in current simulations, we have $\left|\tilde{\mathbf{B}}\right|/B \lt 0.1\%$. In principle, the additional electromagnetic transport can be nonetheless large, depending on the efficiency of parallel (heat) transport along the perturbed magnetic field [54]. However, we have verified that the effect is small in present low-β simulations by additional tests that included electromagnetic flutter (that is therefore otherwise disabled to save a factor 2 in computational time).

Figure 6 shows the difference between one snapshot's density and mean background density – in absolute units of 10¹⁹ m⁻³ (left), as well as normalized to the local mean density (right). Fluctuations are much more prominent on the LFS than on the HFS, which is typical for ballooning modes. In absolute units, fluctuations are negligible in the SOL compared to the confined region – relative to the small local background density, however, SOL fluctuations are much larger, reaching up to 100%, particularly in the region around the X-point where the poloidal magnetic field is weak.

Zoom In Zoom Out Reset image size

While figure 6 clearly demonstrates poloidal asymmetries in the transport, we can still compute average radial diffusivities in quasi-steady state. For E × B particle and heat fluxes

$$\begin{align} \Gamma_r = \left\langle \mathbf{v}_E \cdot \mathbf{e}_r n \right\rangle, \,\,\,\,\, Q^\mathrm{e,i}_r = \left\langle \frac{3}{2} \mathbf{v}_E \cdot \mathbf{e}_r T_\mathrm{e,i} n \right\rangle \end{align} \tag{ 23 }$$

we define diffusivities as suggested by [55] via

$$\begin{align} D_\perp = \frac{\Gamma_r}{\left| \left\langle \partial_r n \right\rangle \right|}, \,\,\,\,\, \chi_\perp^\mathrm{e,i} = \frac{Q^\mathrm{e,i}_r - 1.5 T_\mathrm{e,i} \Gamma_r}{\left\langle n \right\rangle \left| \left\langle \partial_r T_\mathrm{e,i} \right\rangle \right|}. \end{align} \tag{ 24 }$$

In the global average $\left\langle \right\rangle_{r,\theta,\phi,t}$ we obtain for the confined region $D_\perp = 0.18$ m² s⁻¹, $\chi_\perp^\mathrm{e} = 0.27$ m² s⁻¹ and $\chi_\perp^\mathrm{i} = 0.46$ m² s⁻¹ – reasonable values for AUG experiments [22]. $\chi_\perp^\mathrm{i} \gt \chi_\perp^\mathrm{e}$ might imply the presence of the ITG mode – scans of the parallel electron heat conductivity in section 5.2, however, show that this is an effect from diffusive heat flux damping rather than linear instability drive. $\chi_\perp^\mathrm{e} = 2/3 D_\perp$ corroborates the finding that turbulence is driven by ballooning modes [55, 56]. Further evidence could be gained by Fourier analysis in field aligned coordinates, particularly from the parallel envelope and the relative phase shifts between the potential, density and temperatures [51, 52] – but this analysis is complicated for us due to the non-field-aligned grid, and will be deferred to future work.

The presented results also hold in the 1.5 times higher resolution simulations (even higher was not yet computationally feasible, as explained in section 3.1). But it should be remarked, as detailed in the next section, that at higher poloidal resolution the increased zonal flow leads to a radial modulation of profiles and fluctuation levels – a staircase structure [57]. However, in present simulations, this has no impact on overall transport (and profiles, disregarding the modulation) as described in section 3.2.

In this section we want to investigate the composition of the mean radial electric field on closed flux surfaces. Zonal flows will be explained by the dominant contributions to the ion particle balance equation (8), averaged in time and over the flux surface. Additionally, the dominant term in the parallel momentum balance – the ion viscous stress (17) – will reveal a non-trivial contribution from the parallel flow. It is important to note that at any single point in time and space, without averaging, $\partial_t n$ and $\mathbf{v}_E \cdot \nabla n$ are by far the dominant terms in the continuity equation (and similarly for $u_\parallel$ in the parallel momentum equation), but become less important in a large enough ensemble average, consistent with the neoclassical ordering.

As pointed out in the previous section, the electric field is mostly negative in the confined region, as it mainly follows the ion pressure gradient according to $E_r \sim \partial_r p_\mathrm{i}/en$ due to the static equilibrium balance $nC(\varphi) \sim -C(p_\mathrm{i})/e$ (see equation (9)). We are interested in deviations from this balance, and so we show the residue $-C(p_\mathrm{i})/e-nC(\varphi)$ by magenta crosses in figure 8. This residue has to be balanced by other terms in the ion particle balance (8). In the flux surface average, parallel derivatives are annihilated, i.e. $\left\langle \nabla\cdot (nu_\parallel \mathbf{b})\right\rangle_S = 0$, and the E × B advection $\mathbf{v}_E\cdot\nabla n$ is mostly small. Therefore, the residue is mostly balanced by the black curve, which is the divergence of the polarisation particle flux $\nabla\cdot n\mathbf{u}_\mathrm{pol}$ defined in equation (11). In the average, it is a measure of the Reynolds stress. We therefore conclude that the zonal flow can be understood as the residue between E × B and diamagnetic compression which is sustained by the ion polarisation flux (or equivalently the Reynolds stress),

$$\begin{equation*} \left\langle enC(\varphi) + C(p_\mathrm{i}) \right\rangle_{S,t} \approx \left\langle\nabla\cdot(en\mathbf{u}_\mathrm{pol}) \right\rangle_{S,t}. \end{equation*}$$

As explained in the previous section, the zonal flow is driven by near-adiabatic drift waves on Larmor radius scale, and is therefore more pronounced in the higher resolved simulation (but does not necessarily impact overall transport).

Zoom In Zoom Out Reset image size

The remaining deviation between the two curves is explained by hyperviscosity $\mathcal{D}_\Omega$, which is the grid scale numerical dissipation defined in appendix A. It is, however, not just a numerical artifact. As hyperviscosity only acts on grid scale (the value of $\nu_\perp$ is reduced from 3350 to 300 between the 2.5ρ_s0 and the 1.67ρ_s0 resolution simulations), it dissipates the energy arriving there due to the turbulent cascade. This process is enhanced by shear layers, and so dissipation is even stronger in the higher resolution simulation. The separatrix deserves particular attention in this regard, since the E × B shear peaks there, as discussed in the next section. Close to to the separatrix, as marked by the yellow ellipse, this results in a peak of $\nabla\cdot(n\mathbf{u}_\mathrm{pol})$ as well as $\mathcal{D}_\Omega$, and also $\mathbf{v}_E\cdot\nabla n$ largely balanced by $\mathcal{D}_n$.

It is interesting that the divergence of the polarisation particle flux can in fact be conveniently computed from the quasi-neutrality equation (14). In the flux surface average, $\left\langle \nabla\cdot j_\parallel \mathbf{b}\right\rangle_S = 0$, hence $\left\langle \nabla \cdot \left( n \mathbf{u}_\mathrm{pol} \right)\right\rangle_{S,t} = \left\langle C( p_\mathrm{e} + p_\mathrm{i} )/e - \mathcal{D}_\Omega \right\rangle_{S,t}$. This means that the zonal flow results in $\lt C(p)\gt_S \neq 0$, i.e. a (predominantly up-down) pressure asymmetry [60]. Additionally, ballooned transport leads to an inboard-outboard asymmetry [60, 61]. Therefore, neither pressure nor the electrostatic potential can be assumed to be constant along a flux surface ψ (thereby $C(p_\mathrm{e})\approx C(p_\mathrm{i})$).

Let us now discuss the role of ion viscosity. Strictly speaking, the Reynolds stress is only contained in the inertial part of the polarisation particle flux $\nabla\cdot(n\mathbf{u}_\mathrm{pol}^\mathrm{in})$ (see section 2.2). However, we find that the viscous part C(G) is always more than two orders of magnitude smaller than the other contributions in the particle balance, i.e. negligible. This means that in present simulations, viscosity does not directly damp the zonal flow. It does, however, efficiently damp poloidal rotation. To see this, we multiply the parallel momentum balance equation (A3) by density n and the magnetic field strength B and average it in time and over flux surfaces. By this, the parallel pressure gradient is annihilated, $\left\langle B \nabla_\parallel p \right\rangle_S = 0$. The remaining dominant contribution is from the viscous stress function G, or more precisely from its constituents. We have defined the ion viscous stress tensor $\Pi_\mathrm{i}$ in terms of G in section 2.3. We now compare the part of it containing the parallel velocity, plotted in red in figure 8, with the remaining part containing the static equilibrium balance $C(\varphi) + C(p_\mathrm{i})/en$, plotted in blue. As the two parts mostly balance, we conclude that poloidal rotation is near zero according to equation (17),

$$\begin{equation*} \left\langle \mathbf{B} \cdot \nabla \cdot \Pi_\mathrm{i} \right\rangle_{S,t} = 3\eta^\mathrm{i}_0 \left\langle \left( \nabla_\parallel B \right)^2 \right\rangle_{S,t} v_\theta \approx 0. \end{equation*}$$

A small deviation is sustained by $n\mathbf{v}_E \cdot \nabla u_\parallel$ and $p_\mathrm{i}C(u_\parallel)$. Importantly, this condition is fulfilled also in spite of pronounced zonal flows in the high resolution simulation – both parts of the viscous stress tensor are modulated, but cancel each other. This means that while viscosity does not directly damp the zonal flow via C(G), it does damp the resulting mean poloidal rotation by adjusting the parallel velocity $u_\parallel$!

An important observation is that in the lower resolution simulation (2.5ρ_s0), the zonal flow is small, but on average over the radial domain constituents of the stress tensor in the parallel momentum balance are the same. In fact, both simulations develop the same mean toroidal rotation (which is roughly the same as parallel rotation in a large aspect ratio tokamak). Although not strictly valid due to the above mentioned necessary asymmetries along a flux surface, equation (16) suggests that parallel velocity also modifies the electric field locally, which seems to be more important than zonal flows in the low resolution case, figure 7 left, at $\rho_\mathrm{pol}\in(0.94, 0.98)$. The generation, transport and saturation of toroidal rotation are outside the scope of the present study, but a review can be found in [62]. Importantly, in absence of momentum sources in the confined region, as in our present simulations, net toroidal rotation can be only generated in the SOL (see [63]) and transported inward. This is indeed suggested by the parallel velocity profile in figure 5. We point out that ion viscosity seems to be important also for the saturation of this process, as discussed in section 5.1.

The same analysis as above can be performed in a different ensemble average, e.g. toroidal and time (poloidally resolved), as was done in [61] – but more data are required, and results are less concise (2D). A more detailed examination of the generation of mean flows can be found in [41, 64]. We find, consistent with [61], that C(G) is negligible in the vorticity equation. Instead, zonal flows are damped by the GAM oscillation which in turn loses energy through the coupling to Alfvén waves and finally resistivity [48J]. However, we find that this damping is not complete, perhaps because unlike electron pressure, ion pressure does not enter Ohm's law.

We have shown that the electric field is governed by the ion particle balance in the confined region, and by sheath boundary conditions in the SOL – at least in present simulations. The separatrix has its own specific dynamics as it acts as a boundary between closed and open field lines. We can see this by examining what happens to the turbulent eddies. Figure 9 shows the average shearing rate vs. vortex-turn-over rate at the outboard mid-plane. The former is roughly given by $\omega_s = \left| \partial_R \bar{v}_\mathrm{\theta} \right| \approx \left| B^{-1} \partial_R^2 \left\langle \varphi \right\rangle_{t,\phi} \right|$. The latter is estimated by one standard deviation of vorticity $\Omega_\mathrm{std} = \sqrt{\left\langle \Omega_F^2 \right\rangle - \left\langle \Omega_F \right\rangle^2 }$, whereby the fluid vorticity is calculated from the generalised vorticity in GRILLIX via $\Omega_F = \omega_\mathrm{ci0} \hat{B} \hat{\Omega} / \hat{n}$, with $\omega_\mathrm{ci0} = eB_0/M_\mathrm{i}$.

Zoom In Zoom Out Reset image size

Small vortices are elongated, thinned and finally absorbed by the shear flow [4, 58, 65] – driving the zonal flow. Stronger vorticies survive. Importantly, as stated in the last section and visible in figure 8, this process extends up to and peaks at the separatrix. Even though larger zonal flows do lead to somewhat higher shearing, we find that the average vortex-turn-over rate is larger than the shearing rate in the confined region at both 2.5ρ_s0 and 1.67ρ_s0 resolution – explaining why they have barely any effect on overall transport. At the separatrix additional shearing is provided externally due to sheath boundary conditions and the fast SOL outflow, such that the shearing rate exceeds the vortex-turn-over rate and the eddies can be torn apart. In figure 10, at the top of the device, the tilting of eddies is particularly visible. We highlighted one exemplary vortex breaking event (see [66] for a similar experimental finding). Ultimately, vortices are strained out by the shear flow, and dissipated at grid scale. In the particle balance this manifests as a peaked E × B advection rate balanced by dissipation. Supplementary movie 2 shows the dynamics of density fluctuations.

Zoom In Zoom Out Reset image size

Note also how the far SOL structures appear to be nearly frozen in time in contrast – as temperature in hitting ∼3 eV there, time scales are an order of magnitude slower. These structures, pronounced on the scale relative to the background, are in fact tiny due to the very low background $\bar{n} \approx 3\times10^{17}$ m⁻³. Nevertheless, as transport has to saturate globally including the far SOL, its slow dynamics restricts the overall saturation of the simulations.

The viscosity of ions, represented by the viscous stress function G – see sections 2.2–2.3, is an important dissipation mechanism. In section 4.1 we found that C(G), which enters the perpendicular drifts as a higher order correction, is negligible in the polarisation velocity. On the other hand, $\nabla_\parallel G$, which enters the parallel momentum balance at leading order, is crucial for damping poloidal rotation. We have also seen in section 3.3 that this damping leads to significant ion heating.

To closer investigate this, we have conducted a simulation without ion viscosity at all (G = 0). The results at t = 2 ms simulation time are shown in figure 11: outboard mid-plane profiles of density, temperature, parallel velocity and electric field. The density profile is very similar to the reference case in figure 5. The electric field is much larger and almost throughout positive, though. The temperature profiles are broader, fluctuations are larger (up to 15% in density in the confined region) and input power is much larger in this simulation – 482 ± 118 kW for electrons and 525 ± 127 kW for ions, respectively – suggesting increased turbulent transport. The larger fluctuation in input power also suggests stronger GAM oscillations.

Zoom In Zoom Out Reset image size

The electric field is larger in the SOL due to the increased electron temperature gradient. In the confined region, we have performed the same analysis as in section 4.1, finding only a very small contribution from the zonal flow. Although viscosity was disabled in the simulation, we can still compute what the contribution from viscous stress in the parallel momentum balance would have been. Naturally, we find that the viscosity balance (17) is not fulfilled: in the reference simulation, figure 8 left, in an average also over the radial domain, we had $\left\langle \hat{\mathbf{B}} \cdot \nabla \cdot \Pi(u_\parallel) \right\rangle_{S,t,r} = -1.07 \times 10^{26}\, \mathrm{m}^{-2}\,\mathrm{s}^{-2}$ and $\left\langle \hat{\mathbf{B}} \cdot \nabla \cdot \Pi(\varphi,p_\mathrm{i}) \right\rangle_{S,t,r} = 1.10\times 10^{26}$ $\mathrm{m}^{-2}\,\mathrm{s}^{-2}$. In the simulation in this section, without viscosity, one would get $\left\langle \hat{\mathbf{B}} \cdot \nabla \cdot \Pi(u_\parallel) \right\rangle_{S,t,r} = -16.38\times 10^{26}$ $\mathrm{m}^{-2}\mathrm{s}\,^{-2}$ and $\left\langle \hat{\mathbf{B}} \cdot \nabla \cdot \Pi(\varphi,p_\mathrm{i}) \right\rangle_{S,t,r} = 5.23\times10^{26}$ $\mathrm{m}^{-2}\,\mathrm{s}^{-2}$. This means that without viscosity, the plasma obtains a significant poloidal rotation.

At this point, t = 2 ms, the simulation crashes as the electric field starts having large (machine scale) oscillations. A possible reason is the Stringer instability [7, 68], which is expected at low – or absent – viscosity. As the instability is also connected to parallel [60] and toroidal [68] rotation, it is not surprising that we get $\left\langle -\hat{\mathbf{B}} \cdot \nabla \cdot \Pi(u_\parallel) \right\rangle_{S,t,r} \gt \left\langle \hat{\mathbf{B}} \cdot \nabla \cdot \Pi(\varphi,p_\mathrm{i}) \right\rangle_{S,t,r}$. In fact, mean toroidal rotation reaches 26 km s⁻¹ already at t = 2 ms, compared to 13 km s⁻¹ in the saturated reference simulation. However, the inflow of momentum from the SOL [62, 63] and the non-linear Reynolds stress [42] could also be important. At this point, we will not go more into details. But we can conclude that ion viscosity is a crucial mechanism – particularly for damping of poloidal rotation, but also ion heating – and has to be included in realistic tokamak simulations.

The parallel heat conductivity largely determines the parallel heat flux in Braginskii models. Besides obvious consequences for the SOL heat exhaust, circular C-mod simulations [67, 69] have previously shown that flattening of the parallel electron and ion temperature profiles in the confined region can also reduce cross-field transport. To illustrate their impact in current computations, we have repeated the reference simulation with the electron heat conduction lowered to the level of ions, i.e. $\chi_{\parallel \mathrm{e}0} = 35.35$ instead of 940, and sheath heat transmission factor $\gamma_\mathrm{e} = 0$ instead of 2.5. The resulting saturated profiles are displayed in figure 12 (unlike with zero ion viscosity, these simulations were running stably).

Zoom In Zoom Out Reset image size

The most prominent feature is that electron and ion temperature profiles become much broader, including in the SOL, i.e. the pedestal disappears. This is because the profile around the separatrix is determined by the competition of perpendicular transport and parallel outflow, and reducing the heat conductivity and sheath heat transmission significantly hinders the latter. As a secondary effect, due to the increased temperature the ion viscosity increases, damping the parallel flow more strongly. Further, the electric field well deepens. Again, a similar analysis as in section 4.1 reveals that, in this case, this is due to the increased polarisation flux, i.e. zonal flow. This can be either due to a change in the linear instability drive [56], or due to an increased effective resolution: even though the nominal resolution is still 2.5ρ_s0, the separatrix temperature has increased by roughly a factor 2, such that the effective resolution in terms of the local Larmor radius has increased by $\sqrt{2}$. Similarly as with reduced ion viscosity, the density fluctuation level increases up to ∼15% in the confined region, and input power increases to 535 ± 232 kW for electrons and 463 ± 202 kW for ions, respectively. Noticeably, not only the fluctuations increase, but also electrons are now transporting more heat than ions. This corroborates the hypothesis that transport is driven by ballooning rather than ITG modes, except that for the reference case (with high electron heat conduction) the electron heat flux is more suppressed than the ion heat flux by parallel conductivity.