# Generation of zonal flows and impact on transport in competing drift waves and interchange turbulence

Olivier Panico, Yanick Sarazin, Pascale Hennequin, Ozgur Gürcan, Guilhe

Dif-Pradalier, Xavier Garbet, Robin Varennes

► To cite this version:

Olivier Panico, Yanick Sarazin, Pascale Hennequin, Ozgur Gürcan, Guilhe Dif-Pradalier, et al. Generation of zonal flows and impact on transport in competing drift waves and interchange turbulence. Journal of Plasma Physics, 2025, 91 (4), pp.E118. 10.1017/s0022377825100603. hal-05216730

HAL Id: hal-05216730

https://hal.science/hal-05216730v1

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Generation of zonal flows and impact on
transport in competing drift waves and
interchange turbulence

Olivier Panico¹‚²✉, Yanick Sarazin²✉, Pascale Hennequin¹✉, Ozgur Gürcan¹, Guilhe Dif-Pradalier²✉, Xavier Garbet²‚³✉ and Robin Varennes³✉

¹LPP, CNRS, Ecole Polytechnique, Sorbonne Université, Institut Polytechnique de Paris,
Palaiseau 91128, France

²IRFM, CEA, Saint-Paul-lez-Durance F-13108, France

³Nanyang Technological University, Singapore 637371, Singapore

**Corresponding author:** Olivier Panico, olivier.panico@lpp.polytechnique.fr

(Received 28 February 2025; revision received 11 June 2025; accepted 12 June 2025)

The generation and radial structure of zonal flows are studied in competing collisional drift waves and interchange turbulence using the reduced flux-driven nonlinear model TokamID. Zonal flows are generated in both the interchange dominated and adiabatic regimes with the former favoring radially structured flows and avalanche transport. The distance to the instability threshold proves to be key, with a more stable radial flow structure emerging near the threshold and increased energy stored in the flows for interchange turbulence. The avalanches are shown to perturb zonal flow structures in drift-wave turbulence and to reactivate them in the interchange regime. Finally, the ExB staircases with radially structured, stable in time zonal flows are proved beneficial for the overall confinement.

**Key words:** plasma flows, plasma instabilities, plasma simulation

**1. Introduction**

Cross-field transport is critically important for fusion devices such as tokamaks, as it governs the overall energy and particle confinement time. This transport is mainly driven by turbulence at the ion scale. The question of turbulence saturation is then crucial for the design of future machines. The underlying saturation mechanisms also represent a fundamental question, still largely debated. Three main mechanisms can be identified for turbulence saturation in tokamak core plasmas: the relaxation of the profiles (density *n*, temperature *T*), the energy transfer towards dissipative scales and the generation of secondary structures such as zonal flows (ZFs). The second and last mechanisms are partly related since ZFs can contribute to the generation of small-scale structures via vortex shearing.

The profile relaxation is a direct result of the cross-field transport induced by turbulence. A lower thermodynamic gradient reduces the drive for the instability
and resulting turbulence. The transfer towards dissipative scales can occur through mode–mode couplings which lead to turbulent cascades such as those described by Kolmogorov (1941): the energy flows down from the injected scales towards the small scales where it gets dissipated through viscosity. Additionally, the magnetic and velocity shears can lead to the decorrelation of turbulent structures into smaller cells (Biglari, Diamond & Terry 1990; Terry 2000). Finally, the generation of ZFs – large-scale symmetric flows, constant on magnetic surfaces – is considered crucial in tokamaks because their action on turbulence is twofold. First, they are generated by the turbulence and store a part of its energy without inducing cross-field transport. Second, ZFs induce a velocity shear that participates in the shear decorrelation mechanism (Lin et al. 1998). Third, they can also be radially structured, resulting in a set of micro-barriers which is expected to efficiently limit the extension of the radial transport. This process can generate ExB staircases in the form of localized velocity shear layers associated with steps in the pressure profile (Dif-Pradalier et al. 2010, 2015). The transport in between the staircase steps is often considered mediated by avalanches: quasi-ballistic transport events of heat or particles (Dif-Pradalier et al. 2017).

Zonal flows are generated by turbulence through the action of the Reynolds stress – coupling of velocity fluctuations in the radial x and poloidal y directions – being the sum of two contributions, the electric $\Pi_E$ (Diamond & Kim 1991) and the diamagnetic $\Pi_*$ (Smolyakov, Diamond & Medvedev 2000; Hallatschek 2004; McDevitt et al. 2010; Sarazin et al. 2021; Dif-Pradalier et al. 2022): $\Pi_{tot} = \Pi_E + \Pi_*$.

The first, $\Pi_E = \langle \tilde{v}_{EY} \tilde{v}_{Ex} \rangle$, relates to the electric drift $\mathbf{v}_E = \mathbf{E} \times \mathbf{B}/B^2$ which involves the electric $\mathbf{E}$ and magnetic $\mathbf{B}$ fields. The second, $\Pi_* = \langle \tilde{v}_{EY} \tilde{v}_{*x} \rangle$, stems from the diamagnetic drift, $\mathbf{v}_* = \mathbf{B} \times \nabla p/(enB^2)$ which involves the plasma pressure $p$, electron charge $e$ and plasma density $n$.

The generation of ZFs and their impact on turbulence have been studied extensively for the past 20 years (Diamond et al. 2005). While their important role has been assessed, the turbulence regimes leading to their generation and radial structure are still actively explored in view of quantitative prediction. In collisional drift-wave turbulence (Hasegawa & Wakatani 1983), a bifurcation is observed from turbulence dominated regimes to flow dominated regimes. On the one hand, gradient-driven models tend to indicate a collapse of the relative ZF energy at large collisionality (Numata, Ball & Dewar 2007). On the other hand, Panico et al. (2025) with a flux-driven model – where the profiles can adapt freely to the turbulence parameters – found a more gradual reduction of the ZF energy when increasing collisionality. The collisional drift-wave (CDW) instability is controlled using the adiabaticity parameter $C = (k_{\parallel} \rho_s)^2 \omega_{ce} / v_{ei}$, being the parallel wavenumber $k_{\parallel}$ normalized to the Larmor radius $\rho_s = \sqrt{m_i T_i / (eB)}$ divided by the electron-ion collision frequency normalized to the electron cyclotron frequency $\omega_{ce} = eB/m_e$, where $m_i$ and $m_e$ are respectively the ion and the electron mass. However, in tokamaks, it appears important to also take into account the interchange instability that originates from the magnetic field inhomogeneity. The interchange instability is controlled through an effective gravity parameter $g = 2\rho_s/R$, ratio of the Larmor radius to the major radius $R$ of the tokamak. This latter instability has been shown to be important in describing avalanche transport in early reduced flux-driven models (Sarazin & Ghendrih 1998; Garbet et al. 1999). Synergy between interchange and CDWs has also proved to be important in the context of gradient-driven simulations (Scott 2005). Exploring the phase space in this context showed that different confinement properties can be found depending on the turbulence regime (Rogers, Drake & Zeiler 1998; Eich et al. 2021).
In this framework, the goal of the present paper is threefold. First, we use the Tokam1D model derived in Panico et al. (2025) to study the generation of ZFs in a flux-driven competing CDW – interchange turbulence. Second, we investigate the radial structure of the flows and in particular the regimes leading to the appearance of staircases. Third, we explore the impact of ZFs and their radial structure on the transport and confinement. These questions will be tackled by performing several scans of the turbulence parameters C and g on particle confinement time scales while still resolving the small turbulence scales.

In § 2, the Tokam1D model is briefly recalled and the performed simulations are detailed. In particular, the role of the source is discussed. In § 3, the generation of ZFs is studied. The regimes leading to a large ratio of flow to turbulence energy are understood through a study of the underlying Reynolds stress components. Additionally, the ZFs radial structure and the impact of the profile’s distance to the instability threshold are studied. Finally, in § 4, the impact of the flows – and more particularly of their radial structure – on the overall confinement is investigated. Interestingly, depending on the turbulence regime, the ZFs are shown to interact with passing avalanches either constructively or destructively.

## 2. Tokam1D model and performed simulations

Tokam1D is a flux-driven reduced code derived in Panico (2024) and Panico et al. (2025) that features two dominant instabilities of the plasma edge: CDWs and interchange. Assuming an electrostatic isothermal plasma with hot ions, the role of the diamagnetic component of the Reynolds stress is considered. The model stems from the electron density conservation equation and charge conservation equation solved in a slab geometry. It evolves the logarithm of the density $N = \ln n$ and the general vorticity $\Omega = \nabla_{\perp}^2(\phi + \tau N)$ self-consistently, with $\phi$ the electric potential and $\tau = T_i/T_e$ the ion to electron temperature ratio

$$ \partial_t N + \{\phi, N\} = g\partial_y(\phi - N) + \sigma_0 \nabla_{\parallel}^2(N - \phi) + D\nabla_{\perp}^2 N + S_N, \quad (2.1) $$

$$ \partial_t \Omega + \nabla_{\perp,i} \{\phi, \nabla_{\perp,i}(\phi + \tau N)\} = -(1 + \tau)g\partial_y N + \sigma_0 \nabla_{\parallel}^2(N - \phi) + \nu \nabla_{\perp}^2 \Omega. \quad (2.2) $$

Here, $\sigma_0 = \omega_{ce}/\nu_{ei,0}$ is the plasma conductivity taken at a reference plasma density, and D and $\nu$ are dissipation coefficients. Note that the time is normalized to the ion cyclotron frequency $\omega_{cs} = (eB)/m_i$, the lengths to the sound Larmor radius $\rho_s = (m_i c_s)/(eB)$ and the electric potential is normalized to the electron temperature $\phi = e\Phi/T_e$. The model is further reduced to one dimension by separating the quantities into so-called equilibrium (actually flux-surface averages) and fluctuations

$$ N_{eq}(x, t) = \frac{1}{(2\pi)^2} \int_{y,z} N dydz, \quad (2.3) $$

where x stands for the radial, y the poloidal and z the parallel direction. The fluctuating components are Fourier transformed and projected onto single parallel and poloidal wave vectors ($k_{\parallel}$, $k_y$)

$$ \left(\begin{matrix} N \\ \phi \end{matrix}\right) = \left(\begin{matrix} N_{eq} \\ \phi_{eq} \end{matrix}\right) (x, t) + \left(\begin{matrix} N_k \\ \phi_k \end{matrix}\right) (x, t) \exp \left[i(k_y y + k_{\parallel} z)\right] + \text{cc}, \quad (2.4) $$

where cc stands for the complex conjugate. The resulting system solves the equilibrium density $N_{eq}$, the equilibrium velocity $V_{eq} = \partial_x \phi_{eq}$, the density fluctuations $N_k$
and the vorticity fluctuations $\Omega_k$ both taken at the poloidal wavenumber $k_y$. The TokamalD equations read

$$ \partial_t N_{eq} = -\partial_x \Gamma_{turb} + D_0 \partial_x^2 N_{eq} + S_N, \qquad (2.5) $$

$$ \partial_t V_{eq} = -\partial_x \Pi_{tot} + v_0 \partial_x^2 V_{eq} - \mu V_{eq}, \qquad (2.6) $$

$$ \partial_t N_k = +ik_y(\phi_k \partial_x N_{eq} - V_{eq} N_k) + igk_y(\phi_k - N_k) + C(\phi_k - N_k) + D_1 \nabla_{\perp}^2 N_k, \qquad (2.7) $$

$$ \begin{aligned} \partial_t \Omega_k = & -ik_y g(1+\tau)N_k - ik_y V_{eq} \Omega_k + ik_y \partial_x [\phi_k \partial_x (V_{eq} + \tau \partial_x N_{eq})] \\ & - ik_y \partial_x V_{eq} \partial_x (\phi_k + \tau N_k) + C(\phi_k - N_k) + v_1 \nabla_{\perp}^2 \Omega_k. \end{aligned} \qquad (2.8) $$

The reduction to one dimension effectively removes mode-mode (in $k_y$) interaction as those terms involve modes that are not accounted for in the model. The principal consequence is that turbulence saturation through cascade processes is not included in the model. However, nonlinear terms are still retained through the interaction with the equilibrium fields, without any scale separation assumption in time or in the radial direction. As such, turbulence saturation principally occurs through profile relaxation, ZF generation and shear-induced transfer towards stable radial wavenumbers. In this model, it is argued that those mechanisms are essential in fusion plasmas close to marginality (Panico et al. 2025). Additionally, Staebler et al. (2016) have shown that, in the gradient-driven regime, while mode-mode coupling and ZF shearing are competing mechanisms for turbulence saturation at low $k_y$, the latter largely dominates at large $k_y$ (see their figure 1).

The code is flux driven with a particle source $S_N$ and no scale separation in $x$ and time is assumed between equilibrium and fluctuating components. The equilibrium density (2.5) results from the balance between the turbulent flux of particles $\Gamma_{turb} = -2k_y \mathrm{Im}(N_k \phi_k^*)$, where $\mathrm{Im}$ stands for the imaginary part, and the particle source. The equilibrium velocity evolved in (2.6) is that of the ZFs generated through the action of the total Reynolds stress $\Pi_{tot}$ defined as the sum of the electric and diamagnetic components

$$ \Pi_{tot} = \Pi_E + \Pi_* = -2k_y \mathrm{Im}[(\phi_k^* + \tau N_k^*) \partial_x \phi_k]. \qquad (2.9) $$

Each instability is controlled through a physical parameter. The interchange instability is driven by the magnetic curvature $g = 2\rho_s/R$ with $\rho_s$ the sound Larmor radius and $R$ the tokamak major radius. The CDW instability is governed by the adiabatic parameter $C = (k_\parallel \rho_s)^2 \omega_{ce}/v_{ei}$, with $v_{ei}$ the electron-ion collision frequency normalized to the electron cyclotron frequency $\omega_{ce}$ (Hasegawa & Wakatani 1983). Dissipation and viscosity are accounted for via the coefficients ($D_0$, $v_0$) for the equilibrium and ($D_1$, $v_1$) for the fluctuations. An additional friction $\mu$ is considered for the flows, representative of magnetic pumping or neoclassical friction.

The order of magnitude of $C$ and $g$ can be estimated from plasma parameters. Typically, the parallel wavenumber is estimated considering a connection length, $L_q = 2\pi q R$, with $q$ the safety factor, which gives $k_\parallel = 2\pi/L_q = 1/(qR)$. Note that $k_\parallel$ can be characterized experimentally using long range correlation, see for example Mahdizadeh et al. (2007). Estimations of the main parameters are given for standard edge values of the WEST, TCV and MAST-U tokamaks in table 1. The values are computed considering a major radius $R = 2.5$ m and minor radius $a = 0.5$ m for WEST, ($R, a) = (0.87, 0.25)$ m for TCV and ($R, a) = (0.9, 0.6)$ m for MAST-U.

It appears from the parameter dependencies that the strong variation of density and temperature profiles at the edge of tokamak plasmas translates into a wide range
<table><caption>TABLE 1. Main model parameters and their range for typical values of WEST, TCV and MAST-U tokamaks. Plasma parameters are computed at the separatrix: $T_e = T_{sep}$, $n_0 = n_{sep}$ and $B = B_{sep}$. The parallel wavenumber is computed assuming a connection length $k_\parallel = 1/(q_{95}R)$.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Model parameters</th><th colspan="4">Plasma parameters</th></tr><tr><th>g = 2ρ<sub>s</sub>/R</th><th>C = (k<sub>||</sub>ρ<sub>s</sub>)<sup>2</sup>σ<sub>0</sub></th><th>n<sub>sep</sub> [10<sup>19</sup>]</th><th>T<sub>sep</sub> [eV]</th><th>B<sub>sep</sub> [T]</th><th>q<sub>95</sub></th></tr></thead><tbody><tr><td>Dependencies</td><td>α T<sub>e</sub><sup>1/2</sup>/(RB)</td><td>α k<sub>||</sub><sup>2</sup>T<sub>e</sub><sup>5/2</sup>/(n<sub>0</sub>B)</td><td></td><td></td><td></td><td></td></tr><tr><td>WEST</td><td>2–3 × 10<sup>-4</sup></td><td>1–10 × 10<sup>-4</sup></td><td>1–2</td><td>50–100</td><td>3</td><td>5.5</td></tr><tr><td>TCV</td><td>8–12 × 10<sup>-4</sup></td><td>1.5–15 × 10<sup>-3</sup></td><td>0.5–1</td><td>25–50</td><td>1.4</td><td>3.5</td></tr><tr><td>MAST-U</td><td>3 × 10<sup>-3</sup></td><td>3 × 10<sup>-3</sup></td><td>1.5</td><td>40</td><td>0.43</td><td>6</td></tr></tbody></table>

of g and C values. It is interesting to note that C decreases when density increases. It has implication for the interpretation of the nonlinear simulations: a high density (respectively low C) leads to a decrease in the ZF energy for CDW turbulence, see § 3.1.

In order to characterize the generation of ZFs and their impact on the transport as a function of the parameters C and g, a total of 104 simulations are performed starting from a stable density profile with constant dissipation coefficients $D_0 = v_0 = D_1 = v_1 = 10^{-2}$, friction $\mu = 10^{-4}$, poloidal wavenumber $k_y = 0.3$ and ion to electron temperature ratio $\tau = 1$. Simulations are performed on a radial domain of $L_x = 400$ $\rho_s$ using a grid of 1024 points. Dirichlet boundary conditions ensure vanishing fluctuations and vanishing equilibrium velocity on both boundaries. A Neumann boundary condition is used for the equilibrium density on the left-hand side of the simulation and Dirichlet ($N_{eq} = 0.1$) on the right-hand side. Each simulation is run until the density profile has reached statistical steady state. The particle confinement time is defined as

$$ \tau_p = \frac{\int N_{eq} \mathrm{d}x}{\int S_N \mathrm{d}x}. \qquad (2.10) $$

The source is taken Gaussian with the maximum located at the inner boundary of the domain. In order to develop turbulence, one needs to ensure that the maximum gradient achievable by the simulation<fnref n="1" /> – when the diffusion balances the source in (2.5) – called the diffusive gradient $|\partial_x N_{eq}|^{\text{diff}} = \int S_N/D_0$, is larger than the linear threshold defined by the critical gradient $|\partial_x N_{eq}|^{\text{crit}}$. The critical gradient is computed by solving the linear dispersion relation of the Tokam1D system with $V_{eq}^n = 0$ (the equilibrium velocity and all of its derivatives) and $N_{eq}^p = 0$ for $p \ge 2$. Details of the computation are provided in Panico (2024). Two sets of simulations are performed with different source strategies. The first is performed with a constant source $S_N(0) = 10^{-4}$, chosen such that any simulation is above its linear threshold. Since the explored parameter space is large and some cases exhibit large $|\partial_x N_{eq}|^{\text{crit}}$, some simulations yield steady-state gradients much larger than their linear threshold. To palliate this problem and take into account the linear instability threshold depending

<footnote n="1">¹In principle, a local corrugation – in time and space – of the density profile can create a gradient that exceeds the linear instability threshold, triggering turbulence. However, in practice, we rely on the averaged box-size gradient to establish a general threshold.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Diffusive (diamond), critical (square) and steady-state (circle) gradients. The first corresponds to the maximum gradient achievable by the simulations, the second to the linear instability threshold and the last to the mean gradient at steady state. (a) Simulations with a constant $S_N(0) = 10^{-4}$. (b) Simulations with adapted sources such that $|\partial_x N_{eq}|^{\text{diff}} = 6|\partial_x N_{eq}|^{\text{crit}}$. Both figures as a function of C for $g = 10^{-4}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Scanned parameter space for constant source cases, each point corresponds to a simulation that has reached statistical steady state. The color indicates the absolute value of the gradient at the steady state. (circles) scan of C at $g = 10^{-4}$; (triangles) scan of C at $g = 10^{-3}$; (squares) scan of C at $g = 5 \times 10^{-3}$. (diamond) scan of g at $C = 10^{-3}$.</figcaption></figure>

on the turbulence parameters (C, g), a second set of simulation is performed using an adapted source similarly to Panico et al. (2025). In this latter case, the source is chosen such that the diffusive gradient is 6 times the critical gradient. This ensures that each simulation remains relatively close to its linear threshold. Note that the choice $|\partial_x N_{eq}|^{\text{diff}} = 6|\partial_x N_{eq}|^{\text{crit}}$ is arbitrary and has been made such that the turbulent flux $\Gamma_{\text{turb}}$ is larger than the diffusive flux $\Gamma_{\text{diff}}$ at steady state. The diffusive, critical and steady-state gradients are shown in figure 1 for a scan of the parameter C at $g = 10^{-4}$. The case with a constant source is shown on the left, the one with adapted sources on the right.

Performing simulations using a constant source is closer to experiments. However, the effective applied forcing is not constant as a function of the turbulence parameters. Indeed, simulations that have a lower instability threshold – at low C – are forced more than simulations at large C. As a result, the steady-state gradient can be forced very far from its linear threshold, especially for the low C cases. Using an adapted source allows one to have a constant distance to the threshold. Therefore, the system’s stiffness, i.e. the relation between the turbulent flux and the equilibrium density gradient, can be studied. In the following, both results from the constant and adapted sources will be compared.

The scanned parameter space is displayed in figure 2. Each point represents a simulation that has reached statistical steady state. The color indicates the
absolute steady-state gradient. Here, C is scanned from $C = 2 \times 10^{-4}$ to $8 \times 10^{-2}$ for $g = 10^{-4}$, $10^{-3}$ and $5 \times 10^{-3}$ and g is scanned from $10^{-4}$ to $3 \times 10^{-2}$ for $C = 10^{-3}$. The figure features a total of 52 simulations performed at constant source. Large gradients at large C result from the stabilization of both instabilities when getting close to the adiabatic regime. As a matter of fact, simulations at large $C > 6 \times 10^{-2}$ correspond to low turbulence regimes where the turbulent flux is comparable to the diffusive flux. The lowest gradients are obtained for the cases at large g, small C, when the interchange instability dominates. Note that very large values of $g > 10^{-2}$ are stabilizing because of the compressibility terms – divergence of the electric drift and diamagnetic flux in the density conservation equation – and that very low values of $(C, g) = (2 \times 10^{-4}, 10^{-4})$ lead to the stabilization of the CDW instability. Three regimes are identified: the CDW, the interchange and the adiabatic regimes. The transition from one to the other is loosely defined. The CDW dominates at low g medium C, interchange at low C large g and the adiabatic regime takes over at very large values of C regardless of g.

### 3. Generation of zonal flows in competing turbulence

In this section, we study in detail the generation of ZFs in competing CDW-interchange turbulence. The ZFs act on turbulence in two principal ways. The first is by storing energy: the more energy is stored by the ZFs, the less is available for turbulence to produce transport. The second is by inducing a shear that will decorrelate turbulent eddies. The latter mechanism requires ZFs to live on a longer time scale than the typical lifetime of turbulent structures (Biglari et al. 1990; Hahm et al. 1999).

#### 3.1. Flow dominated regimes at large C or g

Numata et al. (2007) have shown that there was a collapse of relative energy stored into the flows at low C. Since $C \propto 1/v_{ei}$, large collisionality is expected to lead to a strong decrease of ZF activity. Later on, Guillon & Gürcan (2025) brought forward a hysteresis mechanism around the transition point. In another contribution, by using a flux-driven code, Panico et al. (2025) have hinted that the decrease of ZF energy could be more gradual than previously expected. Gradient-driven and flux-driven contributions are in agreement provided that one accounts for the evolution of the density gradient together with the adiabatic parameter. However, in all the above cases the turbulence is driven by CDWs only. In this section, we propose to come back to the turbulence–flow energy partition by including the interchange instability. Additionally, by comparing simulations with constant and adapted sources, the impact of the distance to the linear threshold is studied.

To compare the different simulations, the flow-to-turbulence energy partition of the system is used. The energies of the flows and turbulence, $E_{veq}$ and $E_{turb}$ (derived in Panico et al. 2025), are recalled

$$E_{veq} = V_{eq}^2, \qquad (3.1)$$

$$ \begin{aligned} E_{turb} = & 2(1 + \tau)|N_k|^2 + 2|\partial_x \phi_k|^2 + 2|\tau \partial_x N_k|^2 + 4\tau \text{Re}(\partial_x \phi_k \partial_x N_k^*) \\ & + 2k_y^2 [|\phi_k|^2 + |\tau N_k|^2 + 2\tau \text{Re}(\phi_k N_k^*)] \\ & - 2k_y \Im [(\phi_k \partial_x \phi_k^*) + \tau N_k \partial_x \phi_k^* + \tau \phi_k \partial_x N_k^* + \tau^2 N_k \partial_x N_k^*]. \end{aligned} \qquad (3.2) $$

The parameters $E_{veq}$ and $E_{turb}$ are further coarse grained on 4 autocorrelation times of the turbulence $\tau_{turb}$. Then, the root-mean-square (r.m.s.) radial profiles are
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Energy partition between turbulence and flows as a function of C for three scans at different values of g. Cases with a constant source $S_N(0) = 10^{-4}$ are indicated with full lines, those with an adapted source are shown with dotted lines.</figcaption></figure>

computed

$$E_{\text{turb}}^{\text{rms}}(x) = \sqrt{\langle E_{\text{turb}}^2(x, t) \rangle_t}, \quad (3.3)$$

with $\langle \dots \rangle_t$ the time average. Finally, to produce a single estimate per simulation, the energy partition ratio is averaged radially $\langle E_{\text{Veq}}^{\text{rms}}/(E_{\text{Veq}}^{\text{rms}} + E_{\text{turb}}^{\text{rms}}) \rangle_x$.

The results are summarized in figure 3 for the cases $g = 10^{-4}$, $g = 10^{-3}$ and $g = 5 \times 10^{-3}$ as a function of C. Simulations with an adapted source ('fixed' distance to linear threshold) are displayed with dotted lines.

For all considered values of g, a flow dominated regime is found at large C, as expected from previous contributions (Numata et al. 2007; Grander, Locker & Kendl 2024; Panico et al. 2025). Note that, in gradient-driven models, the transition point between flow and turbulence dominated regimes depends on the value of the density gradient. In Numata et al. (2007), the collapse occurs at $C = 0.1$ because the gradient is large: $\kappa = \partial_x \ln n_0 = 0.1$. In comparison, the steady-state box-averaged gradients in the present simulations are lower (see figure 2). Therefore, we expect the transition to occur at lower values of C. Adding the interchange instability modifies the trend at low C. The low g case is analogous to cases without interchange: although not collapsing, the ZF energy is subdominant compared with that of turbulence. When g increases, the flows become more important with approximately 30% of stored energy at $g = 5 \times 10^{-3}$. Additionally, simulations with an adapted source – closer to the linear threshold – display a larger flow to turbulence energy ratio.

The reader should be reminded that the friction exerted on the flows $-\mu V_{eq}$ is set constant in those simulations. However, increasing density also increases the friction coefficient. The large amount of energy captured by the flows at large g and small C may be less pronounced if the friction applied to the flows increases with density (low C). Gianakon, Kruger & Hegna (2002) provide a heuristic expression for this neoclassical coefficient. Simulations taking into account the dependence of C and $\mu$ with density is left for future work.

## 3.2. Both components of the Reynolds stress are crucial depending on turbulence regime

The generation of ZFs and of the associated shear is governed by the total Reynolds stress $\Pi_{tot}$, which is the sum of two contributions: electric $\Pi_E$ and diamagnetic $\Pi_*$.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Radial average of the r.m.s. profile of $\Pi_{tot}$ as a function of $C$ for three values of $g$. (left) Simulations with constant source $S_N(0) = 10^{-4}$. (right) Simulations with adapted sources.</figcaption></figure>

(Diamond & Kim 1991), the latter has been less studied (Smolyakov et al. 2000; Hallatschek 2004; Ivanov et al. 2020; Sarazin et al. 2021; Dif-Pradalier et al. 2022). The parameter $\Pi_{tot}$, being the sum of two contributions, depends on their relative amplitude and whether they are in phase. In this section, both components of the Reynolds stress are studied with the objective to understand the two regimes leading to a large flow to turbulence energy ratio, namely at high $C$ and high $g$.

First, let us take a look at the general behavior of $\Pi_{tot}$ through the correlation and relative amplitude of its two contributions. In figure 4, the r.m.s. value of the total Reynolds stress is shown as a function of $C$ for three different values of $g$ for cases with a constant source (left) and with adapted sources (right). A general trend observed in this figure is that the total Reynolds stress increases with $C$ in most of the parameter domain. Furthermore, at low $C$, the interchange dominated case stands out with a somewhat large Reynolds stress, especially for cases with constant sources, where it decreases until the CDW instability takes over at $C \approx 4 \times 10^{-3}$. The simulations siting closer to their linear threshold, displayed in figure 4(b), depart from the results with constant source in figure 4(a). Indeed, the total Reynolds stress is smaller for the same values of $C$ and $g$ (at the exception of very large $C \le 4 \times 10^{-4}$). That is understood from the linear analysis point of view: as the gradient increases, so does the growth rate, which leads to a larger fluctuation amplitude and resulting Reynolds stress. Interestingly, this behavior is not recovered in the flow-turbulence energy partition in figure 3. This indicates that the turbulence intensity keeps increasing with the gradient, thus leading to a larger $\Pi_{tot}$, but that the generated flows are not able to capture all the added energy. Note that this does not imply that the ZF energy decreases at large gradients: it keeps increasing with the drive $\Pi_{tot}$ but less than the background turbulence (see figure 7).

To get more insight into $\Pi_{tot}$, the correlation between the electric and diamagnetic contributions to the total Reynolds stress is shown in figure 5, together with their relative amplitude. To take into account the possible radial structure of the Reynolds stress, the radial average is computed after having performed the ratio of the r.m.s. profiles.

The three different regimes appear clearly in figure 5. The adiabatic regime, at large $C$, displays the two contributions correlated and in phase, independently of $g$. The interchange regime, at low $C$, shows the contributions in phase opposition, all the more so when $g$ is large. The CDW regime, at medium $C \gtrsim g$, indicates that the diamagnetic contribution is dominant. The amplitude ratio reaches its maximum at a value that depends on both ($C, g$) and the distance to the linear threshold. It then
<figure><img src="image_6.png" /><figcaption>FIGURE 5. (Top) Correlation between electric and diamagnetic contributions to the Reynolds stress $C(\Pi_{\star}, \Pi_{E})$. (Bottom) Amplitude ratio $\langle \Pi_{\star}^{rms} / \Pi_{E}^{rms} \rangle_{x}$. (Left) Simulations with constant source $S_{N}(0) = 10^{-4}$. (right) Simulations with adapted sources.</figcaption></figure>

decreases towards one at large $C$. At low $C$, for cases at medium and large $g$, the electric contribution is dominant, especially as the simulation is close to the linear threshold (figure 5b).

The behavior of the tensors is consistent with the underlying turbulence. On the one hand, getting closer to the adiabatic regime at very large $C$ leads to $N_{k} \sim \phi_{k}$, leading to in-phase and roughly equal Reynolds stress contributions. On the other hand, at low $C$ and large $g$, interchange turbulence is characterized by amplitudes of electric potential fluctuations greater than density fluctuations. This directly impacts the resulting tensors with $|\Pi_{E}| > |\Pi_{\star}|$. Similarly, CDW turbulence in Tokam1D is characterized by $|N_{k}| > |\phi_{k}|$, which leads to a dominant diamagnetic Reynolds stress contribution. The transition towards the adiabatic regime is obtained for lower values of $C$ in cases with adapted sources as compared with fixed sources cases. This is consistent with Gürcan (2024), that puts forward $C/\kappa, \kappa$ being the logarithm of the density gradient, as the correct parameter to study the transition towards adiabaticity. Since simulations with an adapted source tend to have a lower steady-state density gradient, they have a larger $C/\kappa$ for the same values of $C$. Interestingly, the case $g = 5 \times 10^{-3}$ with an adapted source shifts from interchange dominated to adiabatic without transitioning to the CDW dominated regime. The resulting behavior of $\Pi_{tot}$ as a function of $C$ and $g$ can be summarized as follows:

(i) Interchange at small $C$, large $g$: $\Pi_{E}^{rms} > \Pi_{\star}^{rms}$, contributions in phase opposition.

(ii) CDW at medium $C$ (or very low $g$): $\Pi_{E}^{rms} > \Pi_{\star}^{rms}$, contributions weakly correlated. The range of $C$ values for this intermediate regime depends on the forcing (constant vs adapted source), hence on the distance to the threshold.

(iii) Adiabatic at very large $C$: $\Pi_{E}^{rms} \sim \Pi_{\star}^{rms}$, contributions in phase.

One can get further insight into the origin of the correlation between $\Pi_{E}$ and $\Pi_{\star}$ by their analytical study. Bearing in mind that fluctuations can be decomposed into amplitude and phase, $\phi_{k} = |\phi_{k}|e^{i\theta_{k}^{\phi}}$, the electric component of the Reynolds stress reads

$$ \Pi_E = -2k_y |\phi_k|^2 \partial_x \theta_k^\phi. \qquad (3.4) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Real (a) and imaginary (b) parts of ratio $\tau N_k/\phi_k$ for three scans as a function of C. Simulations performed with a constant source $S_N(0) = 10^{-4}$.</figcaption></figure>

Similarly, the diamagnetic component reads

$$ \Pi_{\star} = -2k_y \tau |N_k||\phi_k| \left[ \partial_x \theta_k^{\phi} \cos \Delta\theta_k + \frac{\partial_x |\phi_k|}{|\phi_k|} \sin \Delta\theta_k \right], \quad (3.5) $$

$$ = -2k_y |\phi_k|^2 \partial_x \theta_k^{\phi} \tau \frac{|N_k|}{|\phi_k|} \cos \Delta\theta_k - 2k_y |N_k||\phi_k| \sin \Delta\theta_k \tau \frac{\partial_x |\phi_k|}{|\phi_k|}, \quad (3.6) $$

with $\Delta\theta_k = \theta_k^{\phi} - \theta_k^N$ the cross-phase between the density and electric potential fluctuations taken at the poloidal wavenumber $k_y$. The first term on the right-hand side is proportional to the electrostatic component $\Pi_E$ of the Reynolds stress while the second relates to the logarithmic gradient of the amplitude of the electric potential fluctuations and the turbulent flux of particles ($\Gamma_{turb} = -2k_y|N_k||\phi_k|\sin\Delta\theta_k$). The degree of correlation between the two components of the total Reynolds stress depends on the relative weight of the second term with respect to the first. If the second term is negligible, the two tensors are well correlated. In this situation, $\cos\Delta\theta_k$ determines the sign of the phase coupling, i.e. whether $\Pi_{\star}$ and $\Pi_E$ are in phase or in phase opposition.

The real and imaginary parts of $\tau N_k/\phi_k$ – which correspond respectively to $|\phi_k|/|N_k| \cos \Delta\theta_k$ and $|\phi_k|/|N_k| \sin \Delta\theta_k$ – allow one to understand qualitatively the correlation between $\Pi_E$ and $\Pi_{\star}$. It can be computed both from the nonlinear simulations and from linear analysis. The latter is obtained by performing the linear analysis of the system using the r.m.s. value of the equilibrium density gradient at steady state and discarding other equilibrium parameters: $V_{eq} = \partial_x V_{eq} = \partial_x^2 N_{eq} = 0$. The real part of $\tau N_k/\phi_k$ is shown in figure 6(a), the imaginary part is displayed in figure 6(b) for the cases with a constant source. In both figures, the linear estimate is plotted with dotted lines.

Overall, for both the real and imaginary parts, the linear estimate is close to the nonlinear one even though the analysis is performed without including flow shear. First, in figure 6(b), the amplitude of the imaginary part of the fluctuation ratio is small at large C and at small C for the case large g. Consistently, the correlation between the tensors is large in these regimes. Second, the real part of the fluctuation ratio is negative at low C whatever the value of g, indicating contributions in phase opposition. It becomes positive and roughly equal to one at large C, in agreement with $\Pi_E \sim \Pi_{\star}$ in figure 5(a).

The value of $\Delta\theta_k$ and the sign of $\cos \Delta\theta_k$ depend on the underlying turbulence regime. In particular, the sign of the frequency – or equivalently, phase velocity – is shown to be crucial in determining the sign of the phase relationship between
$\Pi_E$ and $\Pi_\star$. Consider the simple case of the density conservation equation with no compressibility, source or diffusion. In the framework of TokamID, the equation reduces to

$$ \partial_t N + \{\phi, N\} = C(\phi - N). \qquad (3.7) $$

In the limit of very small $C$, the density behaves as a passive scalar advected by the electric drift. Then, the Fourier transform on the fluctuations leads to

$$ \hat{N}_{k\omega} = \frac{\omega_{*e}}{\omega} \hat{\phi}_{k\omega}, \qquad (3.8) $$

with $\omega_{*e} = -k_y \partial_x N_{eq}$ the electron diamagnetic frequency. In the absence of compressibility terms, whenever $C$ is negligible, the density and electric potential fluctuations are always correlated and either in phase or in phase opposition depending on the sign of the frequency with respect to $\omega_{*e}$. It is generally accepted that electron-driven instabilities tend to have a frequency of the sign of the electron diamagnetic frequency while ion-driven instabilities have a frequency of the sign of the ion diamagnetic frequency. In the present case, the CDW instability is driven by electrons and $\omega_{*e}/\omega > 0$. The interchange instability can be driven either by the ions or the electrons, and one cannot conclude which within the present model only. However, it is commonly observed in TokamID that simulations driven by the interchange instability have a frequency of the sign of the ion diamagnetic frequency.

All in all, the behavior of the Reynolds stress depends on the underlying turbulence regime through the correlation and the relative amplitude of its two contributions: electric and diamagnetic. Let us summarize the three different regimes characterized,

(i) Interchange: turbulence dominated by electric potential fluctuations in the ion diamagnetic direction. Leads to electric and diamagnetic contributions to the Reynolds stress in phase opposition with the former larger than the latter.

(ii) CDW: turbulence dominated by density fluctuations in the electron diamagnetic direction. Electric and diamagnetic contributions non-correlated with the latter dominating the dynamics.

(iii) Adiabatic: fluctuations roughly in phase and of the same amplitude. Leads to in-phase contributions to the Reynolds stress of the same amplitude.

### 3.3. ZF structure into staircases in interchange regimes close to the threshold

The second mechanism on which ZFs act on turbulence is by inducing a velocity shear. The determining factor is not only the amplitude of the flows but whether or not they are organized into well-defined sheared layers (Kosuga et al. 2014; Dif-Pradalier et al. 2017). A strong shear is expected to tilt and elongate turbulent structures, leading to their decorrelation, provided that the shear persists longer than the lifetime of the turbulent eddies (Biglari et al. 1990).

The radial structure of ZFs in TokamID changes with the physics parameters ($C$, $g$) and with the distance to the linear threshold. In figure 7, some examples of steady-state flow profiles, as a function of time $t$ and radius $X$, are displayed for different values of $g$. On the top row, the simulations are performed at constant source $S_N(0) = 10^{-4}$ for low $g = 3 \times 10^{-4}$ and large $g = 10^{-2}$. On the bottom row, simulations have the same physics parameter but are performed with adapted sources.
<figure><img src="image_8.png" /></figure>

FIGURE 7. Examples of equilibrium (zonal) velocity $V_{eq} = -\langle E_r \rangle$ at steady state: (a, c) (C, g) = (10⁻³, 3 × 10⁻⁴); (b, d) (C, g) = (10⁻³, 10⁻²). Cases (a) and (b) are computed with a fixed source $S_N(0) = 10^{-4}$ and (c) and (d) are computed with an adapted source, respectively $S_N(0) = 1.57 × 10^{-5}$ and $S_N(0) = 1.76 × 10^{-5}$.

With constant source amplitude, the flows tend to structure radially at larger values of *g* compared with cases with adapted sources. Furthermore, with adapted sources – closer to the linear threshold – the flow radial structures are observed to be more stable in time. This will be emphasized further when discussing figure 8. When less structured, such as in figure 7(a), the flows exhibit a complex behavior with splitting (*X* = 120) and merging (*X* = 60) events. Those events are correlated with large turbulent flux bursts, as will be analyzed in §4.2. Note that simulations with a constant source lead to flow amplitudes approximately two times larger than for the adapted source cases. As a matter of fact, larger forcing for the constant source cases leads to a larger turbulence intensity and a stronger ZF drive.

In the present version of Tokam1D, the flows originating from the plasma radial force balance equation are discarded and only ZFs are considered (2.6). Efforts to include this background radial electric field in the Tokam1D model have already been reported in Panico et al. (2025). The present study focuses on the interplay between turbulence and self-generated flows. Possible impacts of the equilibrium $E_r$ field on the overall dynamics of the system are discussed in §4.2. Also, the turbulence parameters C and g are set constant radially. Therefore, the resulting flows are mostly sinusoidal around zero. One can use their energy spectral density, $S_V(k_x)$, computed from the radial Fourier transform of the equilibrium velocity to study their spatial structure. The radial Fourier transform $\hat{V}_{eq}(k_x, t)$ is computed for $n = 100$ time samples that are then averaged: $S_V(k_x) = \sum |\hat{V}_{eq}|^2/n$. The spectra are displayed in figure 8 for the g-scan with constant and adapted sources.

Two groups of simulations can be distinguished in figure 8(a). The first, at low $g \in [10^{-4} - 10^{-3}]$, is dominated by CDW instability and exhibits a broad extremum around $k_x \approx 0.015$ with amplitude and decay qualitatively similar. The second, at
<figure><img src="image_9.png" /></figure>

FIGURE 8. Energy spectral density of the equilibrium velocity $S_V$ as a function of $k_x$ for different values of the interchange parameter $g$. Each spectrum is the average of 100 independent spectra. (a) Simulations using constant source $S_N(0) = 10^{-4}$. (b) Simulations using adapted sources.

<figure><img src="image_10.png" /></figure>

FIGURE 9. Equilibrium velocity shear $\partial_x V_{eq}$ and equilibrium density gradient $-\partial_x N_{eq}$. (b) Equilibrium velocity shear and effective diffusivity $D_{eff} = -\Gamma_{turb}/\partial_x N_{eq}$. Both are averaged on $\omega_{cs} \Delta t = 3 \times 10^4$ around time $\omega_{cs} t = 1.985 \times 10^6$ and taken from the case $(C, g) = (10^{-3}, 10^{-2})$ with constant source.

large $g > 10^{-3}$, is driven by interchange and exhibits an overall larger spectral density than the first one. For these simulations, the spectrum amplitudes increase together with $g$. This is in agreement with figure 3: ZFs have a larger energy when $g$ increases. Importantly, a peak appears with its maximum shifting from $k_x \approx 0.02$ to $k_x \approx 0.03$ as $g$ increases. The appearance of a peak in the spectrum indicates a clear structuring of the ZFs for $g > 10^{-3}$. Structures have a size between 30 and 50 $\rho_s$ (slightly thinner at large $g$). Simulation with adapted sources, in figure 8(b), always display a peak in the energy spectral density. In most cases, the flows are structured and stable in time with the exception of very small $g$ cases where some merging and splitting events occur. Similarly, larger overall spectrum amplitudes are found for cases at large $g > 3 \times 10^{-3}$. Overall, a more stable in time flow structure is found for simulations close to the linear threshold.

Simulations presenting ZFs radially structured lead to staircases in the form of steps in the density profile. This is illustrated in figure 9 for case $(C, g) = (10^{-3}, 10^{-2})$. The equilibrium density and velocity are plotted as a function of $X$ for $t = 1.985 \times 10^6$. The effective diffusivity $D_{eff} = -\Gamma_{turb}/\partial_x N_{eq}$ is shown in figure 9(b). To remove small-scale effects, both profiles are coarse grained on $\omega_{cs} t = 3 \times 10^4$. The velocity shear oscillates around zero in an almost symmetric pattern. The equilibrium density gradient is negative as a result of the particle source located at $X=0$. When the velocity pattern is stable in time, the density profile is modulated by the shear so
that large (in absolute value) density gradients are located at velocity shear extrema, independently of the sign of the shear. By contrast, all zeros of flow shear coincide with minima of the absolute value of the density gradient, and therefore with maxima of transport. The combined structure of flow shear and ‘steps’ in the density profile corresponds to the staircase regime described by Hornung *et al.* (2016) and Dif-Pradalier *et al.* (2017). Note that there is a negative correlation between the transport – effective diffusivity $D_{eff}$ – and the absolute value of the density gradient. This departure from a standard local flux vs gradient relation – which implies that the gradient causes the flux – has also been pointed out by Villard *et al.* (2013) for temperature gradient corrugations in global gyrokinetic simulations. Seemingly, micro-barriers of transport are generated in location of large shear. From these results only, we expect that the micro-barriers observed in the presence of staircases lead to a better overall confinement. The role of staircases on confinement is further analyzed in § 4.3.

It is interesting to note that, in reduced models, such as TokamID or modified Hasegawa–Wakatani, the generated ZFs have a sinusoidal pattern (see for example Numata *et al.* 2007; Parker & Krommes 2013; Majda, Qi & Cerfon 2018; Grander *et al.* 2024). Therefore, the distance between the staircase steps is similar to the size of the velocity structures. In more complex models, including flux-driven 5D gyrokinetic codes, the distance between staircase steps is larger than the steps themselves and closer to the avalanche size ≈ 40 $\rho_s$ (Dif-Pradalier *et al.* 2017). A possible explanation is the lack of *q*-profile and mode localization effects in those reduced geometries. Including this, one could expect the turbulence, when dominated by instabilities favoring $k_\| \approx 0$, to localize on rational surfaces and then to generate ZF structures accordingly (Dominski *et al.* 2017). This has been hinted at experimentally on the TJ-II stellarator by Van Milligen *et al.* (2022). Finally, one can mention the model – based on the wave kinetic equation accounting for the dynamics of drift waves coupled to an equation governing the ZF dynamics – proposed by Garbet *et al.* (2021) where different shapes of staircase patterns can be obtained. There, it is found that staircases result 'from the interaction between propagating wave packets (understood as avalanches) and waves that are trapped in zonal flow velocity wells'. The detailed characteristics of the staircases (amplitude, shape and periodicity) are determined by those of the background fluctuations, notably their spectra and growth rates. In particular, structured ZFs appear to exhibit non-sinusoidal radial profiles, peaked at their maxima (located at the O-points of the islands that trap the drift waves in their phase space), when the growth rate of the drift waves is maximum at vanishing radial wavenumber. Experimentally, work from Choi *et al.* (2024) hints towards a Fréchet-like distribution of the staircase steps, in agreement with Dif-Pradalier *et al.* (2017).

To sum up, ZFs are found to be radially structured in near-marginal, interchange-driven plasmas. For all studied cases, a stable and radially localized flow shear always leads to staircase steps in the density profile. The sign of the zonal shear does not play a role and the density steps follow the flow sinusoidal pattern. This is associated with the formation of a micro-barrier of transport in the form of a reduced effective diffusivity.

### 3.4. Larger flows close to the instability threshold

In figure 3, the distance to the linear threshold plays a role in the energy partition of the system. In this section, the objective is to compare two simulations that share
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Relative turbulent flux $\Gamma_{turb}/\Gamma_{tot}$, energy partition $E_{V_{eq}}/(E_{V_{eq}} + E_{turb})$ and distance to linear threshold as a function of the source $S_n$: (a) CDW driven; $g = 10^{-4}$, $C = 2 \times 10^{-2}$ and (b) interchange driven; $g = 5 \times 10^{-3}$, $C = 4 \times 10^{-4}$.</figcaption></figure>

similar flow with turbulence energy ratios but that have different control parameters and hence different turbulence properties. The first, driven by CDW, has a large $C = 2 \times 10^{-2}$ and small $g = 10^{-4}$, the flows are little structured with meandering and splitting events. The second is driven by interchange. It has a large magnetic curvature $g = 5 \times 10^{-3}$ and small adiabaticity parameter $C = 4 \times 10^{-4}$, the flows are highly structured and stable in time.

To vary the distance to the instability threshold, the source is slowly lowered, ensuring that the simulation reaches a steady state at each step. When reducing the source, the simulation approaches its nonlinear threshold, until turbulence is lost. The relative distance to the linear threshold is quantified with the following dimensionless parameter $\Delta_{lin}$:

$$ \Delta_{lin} = \left| \frac{\partial_x N_{eq} - \partial_x N_{eq}^{crit}}{\partial_x N_{eq}^{crit}} \right| \qquad (3.9) $$

To monitor the amount of 'turbulence', the ratio between the turbulent and diffusive fluxes is computed along with the flow-turbulence energy partition. The result is shown in figure 10 as a function of $\Delta_{lin}$. The horizontal error bars account for the corrugation of the equilibrium density gradient, the vertical error bars indicate the standard deviation for each dataset.

As the simulations get closer to the linear threshold, the turbulence intensity decreases and so does $\Gamma_{turb}/\Gamma_{diff}$: turbulence generates a smaller part of the total transport. An important point to notice is the difference in stiffness for drift-wave and interchange turbulence, the latter leading to a larger turbulent flux just above the linear threshold. In both cases, the nonlinear threshold is very close to the linear critical gradient and no significant upshift is found. However, the simulations being flux driven, it is difficult to define properly the gradient and the threshold because the profiles can be corrugated. In those cases, larger gradients can occur locally and produce turbulence even though the averaged profile is below the linear threshold. In agreement with figure 3, the case at (large $g$, small $C$) (figure 10b) displays an increased flow to turbulence energy ratio as the threshold is approached. In this case, the flows store more energy close to the threshold.

The marginality – loosely defined as the distance to a linear/nonlinear instability threshold – has already been considered important for both simulations Dif-Pradalier et al. (2017) and experiments Hornung et al. (2016). In those contributions, it is
reported that staircases are found in near-marginal ion drift-wave turbulence and are not observed in the electron drift-wave regimes. In the present work, the statistics offered by the large number of simulations stands in agreement with those previous observations as the staircases are found in interchange close to the linear threshold.

## 4. Role of zonal flows on transport and avalanches

Ultimately, the generation of ZFs is crucial for a fusion machine because they mitigate the turbulent transport (Lin et al. 1998). Two mechanisms are at play: they constitute a sink of energy for turbulence (and create no transport themselves) and they induce a shearing (Biglari et al. 1990; Diamond & Kim 1991). In particular, the existence of radially localized time stable shear layers of flows is expected to lead to a better confinement (Dif-Pradalier et al. 2017). However, the transport itself is also complex and self-organizes into various structures. In TokamlD, the linear analysis indicates that the growth rate is maximum for the lowest values of $k_x$ accessible by the system. As a result, turbulence tends to develop radially extended vortices, often called streamers, that mix large portions of the plasma. The flow shear is instrumental in breaking those extended structures into smaller cells. Transport also exhibits complex structures in the form of avalanches: almost ballistic transport events of heat and particles on long distances with respect to the turbulence correlation lengths. A single, large transport barrier resulting from a large shear could possibly stop the transport event (Fedorczak et al. 2009). However, the question of the interplay between avalanches and staircases is still debated. Staircases can be expected to either limit the radial extent of avalanches or arise as a result of avalanches themselves, see respectively work from Choi et al. (2024) and Kosuga et al. (2014). In this section, the transport is analyzed in the drift-wave and interchange regimes. It is shown that the confinement largely depends on the control parameters C and g. Then, the role of avalanches is analyzed. They prove to be of importance both as a means to carry particles and as quantities that interact with the dynamical zonal structures. Finally, the influence of staircases on transport is assessed.

First, let us analyze the confinement in a statistical sense by looking at the confinement time of particles (2.10). Note that this quantity yields similar information to the ratio of the turbulent to diffusive fluxes of particles at the steady state. In the framework of TokamlD, the ‘diffusive flux’ is linked to the term $D_0 \partial_x^2 N_{eq}$ in (2.5). At the steady state, a better confinement is obtained when the diffusive flux of particles accounts for a larger part of the transport. Actually, the theoretical ‘optimum’ confinement is when there is no turbulence, hence $\Gamma_{turb} \approx 0$ and $\partial_x N_{eq} \approx \partial_x N_{eq}^{diff}$. The particle confinement time at the steady state is displayed in figure 11.

The particle confinement time $\tau_p$ relates to the steady-state density profile and the source. It is mostly governed by the instability threshold and the system’s stiffness, i.e. its ability to depart from the instability threshold. The threshold – defined loosely as marginality – has already been pointed out as crucial for the overall confinement by Diamond & Hahm (1995). In figure 11, the confinement times of the constant $S_N$ cases are mostly governed by the linear threshold. For the adapted source cases, the linear instability threshold is taken into account in the source, and the confinement depends on the system’s stiffness only. As such, it is easier to infer the role of nonlinear quantities, such as ZFs or avalanches, on the particle confinement time.

In figure 11, the confinement time mostly increases with C for two main reasons: the critical gradient increases with C (cf. figure 1) and the stiffness reduces with C. The threshold is also larger at small C, small g, leading to an increased confinement. Note that the cases with an adapted source yield a larger confinement time
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Confinement time of the particles $\tau_p$ as a function of $C$ for three values of $g$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Examples of turbulent flux of particles $\Gamma_{turb} = -2k_y \tilde{\delta}(N_k \phi_k^*)$ at steady state. (a) $(C, g) = (10^{-3}, 3 \times 10^{-4})$. (b) $(C, g) = (10^{-3}, 10^{-2})$. Both are computed with a fixed source at $S_N(0) = 10^{-4}$.</figcaption></figure>

as they are closer to their linear threshold. They exhibit a larger flow to turbulence energy ratio and a more stable in time radial structuring of ZFs. When the source is larger, the excess of energy is stored into the turbulence which leads to an increased transport.

4.1. Transition to avalanche-like transport in interchange turbulence

The turbulent transport is shown to transit towards avalanche-like transport when increasing $g$ at fixed $C$. Choosing the cases (a) and (b) of figure 7, the turbulent flux of particles as a function of $X$ and time is displayed in figure 12 with a low $g = 3 \times 10^{-4}$ on the left and a large $g = 10^{-2}$ on the right.

At larger $g$, the turbulent flux of particles displays diagonal stripes across large portions of the simulation domain. Those correspond to ballistic transport events of particles compatible with avalanches. Note that the information of the avalanche travels both in and outwards while the flux is always positive, i.e. the transport of particles is always outwards from the domain. The propagation of the avalanche is often understood similarly to a sand pile avalanche (Diamond & Hahm 1995). A local relaxation of the profile leads to two high gradients on each side of the plateau. At those locations, the larger gradient induces a larger growth rate and a resulting increased transport, which leads to a flattening of the steep regions. As such, the perturbation moves in two directions: a bump propagates down gradient while a void
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Propagating avalanches on the equilibrium density profile, case (C, g) = (10⁻³, 10⁻²) with constant source $S_N(0) = 10^{-4}$.</figcaption></figure>

propagates up gradient. Actually, it is possible to verify on the equilibrium density profile whether there is a propagation of voids and bumps by plotting $N_{eq}(x, t) - \langle N_{eq} \rangle_t(x)$. This is done in figure 13 for the case (C, g) = (10⁻³, 10⁻²). The diagonal stripes indicative of avalanches appear clearly. This time, events traveling to the right-hand side are positive perturbations – bumps – in red while events traveling to the left are negative – voids – in blue.

Note that the avalanche patterns exhibit an almost symmetric bi-directional behavior, i.e. as many avalanches propagate inward and outward at each radial position. The direction of propagation of avalanches has been shown to depend on the sign of the mean flow shear in earlier publications. Different mechanisms were proposed to understand this behavior by Idomura et al. (2009) and McMillan et al. (2009). The reason for equi-probable inward- and outward-moving avalanches in TokamID likely stems from the absence of equilibrium radial electric field in the present version of the code – the radial force balance is not enforced. As a result, the radial derivative of $E_r$ does not have a well-defined sign in large radial regions of the system. Conversely, it oscillates with positive and negative values in response to turbulent bursts. Studying the impact of such equilibrium $E_r$ on avalanche dynamics with TokamID is left for future work. The avalanche-like pattern presented above occurs in simulations with highly structured flow profiles. In the observed simulations, when the ZFs are structured, there is always an avalanche-like dominated transport. The avalanches are also visible directly on the flow pattern, as they perturb it when propagating across. More details regarding the interaction with flows are given in the next section.

Note that the two cases presented figure 12 are at both ends of the g-scan. The intermediate values of g lead to intermediate behaviors with avalanches that propagate on smaller parts of the domain. Avalanches can be studied statistically by computing the skewness S and kurtosis K, respectively third and fourth moments of the probability density function (pdf) p of the turbulent flux of particle, $Z = \Gamma_{turb}(t)$

$$S(Z) = E \left[ \left( \frac{Z - \mu}{\sigma} \right)^3 \right] \quad ; \quad K(Z) = E \left[ \left( \frac{Z - \mu}{\sigma} \right)^4 \right], \qquad (4.1)$$
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Probability distribution function of the turbulent flux for simulations at constant (a) and adapted (b) sources as a function of g. All the simulations are computed for C = 10⁻³. Statistics are computed at each radial position.</figcaption></figure>

with $\mu$ the statistical mean and $\sigma$ the standard deviation of the fluctuating sig-
nal $Z(t)$ and $E[F(Z)] = \int_{-\infty}^{\infty} F(Z)p(Z)dZ$. These functions measure respectively
the pdf asymmetry and width (weight of the tails). For a Gaussian statistics, $S=0$
and $K=3$. In the present case, avalanches are radially extended and ‘rare’ events.
Therefore, we expect them to lead to important values for skewness and kurtosis.
For both cases with constant and adapted sources, $S$ vs $K$ is plotted in figure 14 for
the $g$-scans. Simulations with a constant source are indicated in red (left), the ones
with adapted sources are in blue (right).

The kurtosis increases with the skewness in an almost polynomial way. Interestingly, the relation between kurtosis and skewness is similar to the one reported in experiments Labit *et al.* (2007). Also, both the skewness and the kurtosis are shown to increase with larger values of *g*, displaying a more skewed and deviated pdf. Overall, the avalanches account for a larger part of the statistics as *g* is increased and as the system is far from the linear threshold. Note that this behavior is not universal, other contributions for example from Diamond & Hahm (1995) and Sarazin *et al.* (2000) discuss the transport to be avalanche-like close to the linear threshold.

4.2. Interplay between avalanches and zonal flows

As stated in the previous sections, staircases are observed mainly at large *g* where
they coexist with avalanches. To get more insight into the interaction between flows
and avalanches, we plot on top of the flow profiles the 90 % quantile of the turbulent
flux presented in figure 12. In other words, we add the largest events of the turbulent
flux of particles. The result is shown in figure 15 for the case at small *g*, dominated
by CDW and in figure 16 for the case at large *g* dominated by interchange.

In the first case, there are no avalanches *per se* but bursts of turbulent flux of
particles occuring throughout the simulation. The radial extent of those bursts is
roughly of the size of the flow structures themselves. Interestingly, the locations of
the bursts are correlated to modifications of the flow topology. This is the case both
for merging events (bottom left of the zoom) and for splitting events (top right).

On the second case, in figure 16, the avalanches are large enough so that they
propagate across the ZF structures. Here, two types of ZFs can be identified. The
first type consists of ZFs at finite frequency moving along the almost ballistic trajec-
tories of the avalanches of flux: these avalanches leave a faint imprint on the flow
pattern as can be seen on the color map. These ones exhibit the usual predator–prey
<figure><img src="image_16.png" /><figcaption>FIGURE 15. (left) Example of flow as a function of $X$ and time with super-imposed 90% quantile of the particle flux. Case $(C, g) = (10^{-3}, 3 \times 10^{-4})$ with constant source $S_N(0) = 10^{-4}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 16. (left) Example of flow as a function of $X$ and time with super-imposed 90% quantile of the particle flux. Case $(C, g) = (10^{-3}, 10^{-2})$ with constant source $S_N(0) = 10^{-4}$. (right) Reactivation of an existing ZF structure by passing avalanches. Temporal slice taken at $X = 305$.</figcaption></figure>

dynamics (see Diamond *et al.* 1994 for a heuristic model, and e.g. Morel, Gürcan & Berionni 2013 for its observation in gyrokinetic simulations): ZFs feeding on turbulence and contributing to its saturation. This mechanism is suspected to govern the local time duration of an avalanche burst: the flux increases locally, leading to the increase of ZFs with some time delay; the latter then act to reduce turbulence and the associated flux. The second type is made of the structured ZF layers at almost zero frequency which exhibit a peculiar behavior when an avalanche passes through. To investigate the interaction, we display a slice at $X = 305$ as a function of time on the right-hand side of figure 16. The flows decay exponentially in the absence of any turbulent drive. This results from the imposed friction ($-\mu V_{eq}$) in (2.6). Conversely, when an avalanche passes through the layer, the flow structures are observed to be reactivated. The particularity here is to have a staircase structure, where ZFs are radially localized, that is maintained by the passing avalanches. In these localized layers of strong ZFs, the flows (predator) are no longer tied to the avalanche
(prey) dynamics as would be expected in a predator–prey model. Such a regime can be expected when the system is efficient in converting turbulent energy into flow energy (large predator population growth for a given prey population) and/or when the life time of the flows (small damping) is long with respect to the time interval between avalanche bursts (regeneration time of the prey). Note that this kind of dynamics where bursts of turbulence, or proxy for avalanches, contribute to the self-sustainment of long-lived ZFs layers has already been hinted at in both reduced gradient-driven (Ivanov et al. 2020) and gyrokinetic (Dif-Pradalier et al. 2017) models. This generation mechanism has been considered by Zhu, Zhou & Dodin (2020) for a Hasegawa–Wakatani type of model and by Kosuga et al. (2014). In the latter contribution, the key element is to consider a delay between the instantaneous and the mean heat flux, thus leading to ‘jams’ of heat flux waves. A possible way of closing the gap with their model would be to compute the distribution of the avalanches ‘waiting time’ in the present simulations. This comparison is left for future work.

### 4.3. Staircases improve the overall confinement

Do staircases improve the overall confinement? As stated in figure 11, the confinement is mainly governed by the linear threshold and the system’s stiffness. While the effect of the linear threshold is ‘neutralized’ in simulations with the adapted source, the system’s stiffness results from both linear – instability growth rate, linear cross-phase – and nonlinear – ZFs, avalanches – properties. Therefore, it is necessary to normalize the particle confinement time to a linear estimate so that simulations with different turbulence parameters can be compared.

In the framework of the quasilinear theory, the transverse transport is diffusive, characterized by the diffusion coefficient $D_{QL} = \gamma_k \sin(\Delta\varphi_k)/k_y^2$, with $\gamma_k$ the linear growth rate and $\Delta\varphi_k$ the linear cross-phase between density and electric potential fluctuations, both taken at the poloidal wavenumber $k_y$. Therefore, the confinement time is simply defined as the square of the radial distance over which confinement is required, namely $L_x$, divided by this transport coefficient $D_{QL}$. It results that simulations with a large growth rate and cross-phase lead to a small quasilinear confinement time. The chosen estimate reads

$$ \tau_{QL} = \frac{k_y^2 L_x^2}{-\gamma_k \sin \Delta\varphi_k} \qquad (4.2) $$

The minus sign comes from the fact that a negative $\sin \Delta\varphi_k$ leads to an outward turbulent flux $\Gamma_{turb}$. The linear quantities $\gamma_k$ and $\sin \Delta\varphi_k$ are computed without taking into account profile corrugations: only the averaged steady-state equilibrium density profile is considered. However, simulations that include staircases usually also include density profile corrugations. Those lead to larger gradients locally and can give a larger estimation of $\gamma_k$. However, the linear quantities are also computed without taking into account any kind of equilibrium flows: $V_{eq} = \partial_x V_{eq} = 0$, because the objective is to assess the role of the flows on the confinement. All in all, we choose to keep the linear estimate simple by taking only a mean estimate of the gradient and no flows. In the situation where the growth rate or the cross-phase vanishes, the mixing length estimate for the confinement time diverges as no transport is expected. In practice, we choose $\tau_L = \min(\tau_{QL}, \tau_{diff})$ with $\tau_{diff} = \int N_{eq}^{diff} dx / \int S_N dx$ the maximum confinement time achievable in a Tokamalp simulation. The particle confinement time normalized by the linear estimate is shown in figure 17 as a
<figure><img src="image_18.png" /><figcaption>FIGURE 17. Normalized particle confinement time $\tau_p/\tau_L$. (a) As a function of C, note that $\tau_L = \min(\tau_{QL}, \tau_{diff})$. (b) As a function of g for C = $10^{-3}$. Both for simulations with an adapted source $S_N(0) = 10^{-4}$.</figcaption></figure>

function of C and g. The normalization is performed on cases with an adapted source so that the forcing is the same across the parameter space.

In figure 17, the trend is different from the non-normalized particle confinement time $\tau_p$ in figure 11. Simulations leading to the largest $\tau_p/\tau_L$ are those at large values of g. The normalized confinement time also increases with C at low and medium g but then reaches a plateau. We argue that, in the framework of Tokam1D, the normalized confinement time can be used to study the role of ZFs in improving the confinement. To this end, we compare the above figure to the system's energy partition, figure 3, and to the flows radial structure figure 8. The largest flows in terms of energy are found either at large value of C or at large value of g. The second case also leads to radially structured ZFs that induce a more stable shear. The normalized particle confinement time is larger for cases yielding a large flow to turbulence energy ratio and even more so when they are strongly radially structured (large g, small C). This supports the conclusion that flows, and in particular radially structured flows, lead to a better overall confinement.

## 5. Conclusion

Three important topics have been assessed in this paper: (i) the generation of ZFs in competing CDWs – interchange turbulence, (ii) their radial structure and (iii) the role of ZFs and staircases in transport and avalanches.

To this end, the reduced flux-driven nonlinear model Tokam1D that features both CDW and interchange instabilities has been used. The CDW instability is controlled by the adiabatic parameter C which scales like the square of the turbulence parallel wavenumber divided by the electron-ion collision frequency. Interchange is driven by an effective curvature g which, in the framework of Tokam1D, relates to the strength of the magnetic field inhomogeneity. The isothermal and electrostatic model includes a finite ion temperature, so that both the electric and diamagnetic contributions to the Reynolds stress are taken into account in the generation of ZFs. This model is reduced to one dimension by separating the density and electric potential into equilibrium – to be understood as flux-surface averages – and fluctuating components. The latter is Fourier transformed and projected onto single poloidal $k_y$ and parallel $k_\| $ wavenumbers. As such, while mode-mode nonlinear coupling between finite $k_y$ modes is discarded, nonlinear transfer to different radial wavenumbers $k_x$ still occurs through the coupling with the equilibrium fields. Importantly the model is flux driven and does not assume any scale separation between equilibrium and
fluctuating components. This is crucial for the density and flow profiles to corrugate and form staircases (Panico et al. 2025). A total of 104 simulations have been evolved on the particle confinement time while still resolving the turbulence time and spatial scales. Two sets of simulations with different strategies for the source are used with the same parameters C and g. The first is set with a constant amplitude at $S_N(0) = 10^{-4}$ throughout the whole scanned parameter space. The second is adapted for each simulation so that each one remains at a controlled distance from the linear stability threshold.

Three main regimes are identified. The first, at low C and large g, is the interchange regime. The second, at medium C and low g, is the CDW regime. The third, at very large C, is the adiabatic regime. Flows are shown to be generated and to efficiently store a significant part of the turbulence energy in both the interchange and the adiabatic regimes. Unlike previous publications that include only CDW instability, the presence of the interchange instability allows ZFs to still be generated even when the adiabatic parameter is small. One of the consequences is that the ‘ZF collapse’ routinely observed in gradient-driven codes featuring the sole CDW instability, is not present in the TokamD simulations. The interchange dominated regime also exhibits radially structured ZFs – intimately associated with corrugations in the equilibrium density profile – and avalanche transport.

The two regimes leading to large flow-to-turbulence energy ratio are understood from the role of the underlying Reynolds stress. The interchange dominated regime is dominated by the electric contribution to the Reynolds stress with the subdominant diamagnetic contribution in phase opposition. As expected, the adiabatic regime is characterized by similar amplitude and in-phase electric and diamagnetic contributions. The third regime, dominated by CDW instability, is driven by the diamagnetic contribution to the Reynolds stress reaching roughly two times the electric contribution.

The distance from the equilibrium density profile to the instability threshold plays a role in generating the ZFs. First, flows tend to be more structured when close to the instability threshold. Second, in interchange-driven plasmas, flows are able to capture more turbulent energy close to the threshold. This second characteristic is not observed in CDW-driven turbulence for the explored parameter range.

Finally, the role of ZFs in transport is assessed. First of all, the zonal structures are shown to interact with the avalanches. In the CDW dominated case, the turbulent flux of particles is shown to disturb the ZF structures, leading to modifications of their topology. In the interchange-driven case, the large avalanches are shown to reactivate the radially localized ZF structures – otherwise exhibiting an e-folding decay in time due to collisional friction – by traveling through them. At last, using a normalization that takes into account the linear properties of the system, the cases with radially localized ZFs – in the form of staircases – exhibit a larger normalized particle confinement time. Staircases are shown beneficial for the overall confinement.

The results and trends reported in the present contribution, including the generation of ZFs and staircases in near-marginal interchange dominated turbulence, should be understood as guides to conduct simulations and analyze the results of more advanced models, which can ultimately validate these findings. The minimal requirement for these models is to be flux driven and to include both contributions – electric and diamagnetic – to the Reynolds stress.

We can think of at least two main physical ingredients that are missing in the present TokamD model – apart from the nonlinear mode-mode coupling already
discussed at length – and that may affect the reported dynamics. The first one deals with geodesic acoustic modes, the high-frequency branch of the ZFs, that are known to exchange energy with both turbulence and low-frequency ZFs (Miki & Diamond 2010). Actually, there may be a way to incorporate their dynamics, at least heuristically, in an extended version of TokamID by adding another mode at a different $k_y$. The second one relates to the intrinsically three-dimensional nature of turbulence in tokamak plasmas. Especially, higher dimensions allow one to explore the correlated role of the ballooning character of interchange turbulence and of the radial localization of modes around rational flux surfaces, hence the role of the safety factor *q*, possibly impacting the flow structures themselves, as recently suggested (Van Milligen et al. 2022).

Regarding the TokamID model itself, obvious improvements can be envisaged to enrich the addressed physics, hence the relevance of the model. Adding the temperature equations would enable more experimentally relevant modes such as resistive ballooning and ion temperature gradient modes. This would also allow one to study the conditions for the possible decoupling between particle and heat transport channels – some works point towards the critical role the parallel electron heat flux (Manz et al. 2020). Also, describing high-confinement regimes such as the H-mode requires us to take into account electromagnetic effects which are known to affect the electron response in large gradient regimes (Scott 1997; Zholobenko et al. 2024). Also, this would permit us to address the role of additional instabilities such as kinetic ballooning modes. Last, TokamID can also be extended to the scrape-off layer by providing an adequate closure for the parallel dynamics. These improvements are ongoing and will be reported in future works.

## Acknowledgments

The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme 'Anti-diffusive dynamics: from sub-cellular to astrophysical scales', where work on this paper was undertaken.

Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.

## Funding

This work was supported by EPSRC grant EP/R014604/1. G.D.P. and Y.S. acknowledge that this work was partially supported by a grant from the Simons Foundation. This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion). The French contribution to this work has been partially funded by the Research Federation Fusion par Confinement Magnétique (FRFCM). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI. Neither the European Union nor the European Commission nor FRFCM can be held responsible for them. This work was supported by the EUROfusion Theory and Advanced Simulation Coordination (E-TASC) initiative under the TSVV (Theory, Simulation, Verification, and Validation) 'Physics of the L-H Transition and Pedestals' (2021–2025).
**Declaration of interests**

The authors report no conflict of interest.

**Data availability**

The data that support the findings of this study are available upon reasonable request.

**REFERENCES**

BIGLARI, H., DIAMOND, P.H. & TERRY, P.W. 1990 Influence of sheared poloidal rotation on edge turbulence. *Phys. Fluids B* **2**, 1–4.

CHOI, M.J., et al. 2024 Mesoscopic transport in kstar plasmas: avalanches and the e × b staircase. *Plasma Phys. Control. Fusion* **66**, 065013,

DIAMOND, P.H. & KIM, Y.-B. 1991 Theory of mean poloidal flow generation by turbulence. *Phys. Fluids B* **3**, 1626–1633.

DIAMOND, P.H., LIANG, Y.-M., CARRERAS, B.A. & TERRY, P.W. 1994 Self-regulating shear flow turbulence: a paradigm for the l to h transition. *Phys. Rev. Lett.* **72**, 2565–2568.

DIAMOND, P.H. & HAHM, T.S. 1995 On the dynamics of turbulent transport near marginal stability. *Phys. Plasmas* **2**, 3640–3649.

DIAMOND, P.H., ITOH, S.-I., ITOH, K. & HAHM, T.S. 2005 Zonal flows in plasma – a review. *Plasma Phys. Control. Fusion* **47**, R35–R161.

DIF-PRADALIER, G., DIAMOND, P.H., GRANDGIRARD, V., SARAZIN, Y., ABITEBOUL, J., GARBET, X., GHENDRIH, P., STRUGAREK, A., KU, S. & CHANG, C.S. 2010 On the validity of the local diffusive paradigm in turbulent plasma transport. *Phys. Rev. E* **82**, 025401.

DIF-PRADALIER, G., et al. Transport barrier onset and edge turbulence shortfall in fusion plasmas, *Commun. Phys.* **5**, 229.

DIF-PRADALIER, G., HORNUNG, G., GARBET, X., GHENDRIH, PH, GRANDGIRARD, V., LATU, G. & SARAZIN, Y. 2017 The e × b staircase of magnetised plasmas. *Nucl. Fusion* **57**, 066026.

DIF-PRADALIER, G., et al. 2015 Finding the elusive e × b staircase in magnetized plasmas. *Phys. Rev. Lett.* **114**, 085004.

DOMINSKI, J., MCMILLAN, B.F., BRUNNER, S., MERLO, G., TRAN, T.-M. & VILLARD, L. 2017 An arbitrary wavelength solver for global gyrokinetic simulations. application to the study of fine radial structures on microturbulence due to non-adiabatic passing electron dynamics. *Phys. Plasmas* **24**, 022308.

EICH, T., et al. 2021 The separatrix operational space of ASDEX Upgrade due to interchange-drift-Alfvén turbulence. *Nucl. Fusion* **61**, 086017.

FEDORCZAK, N., GUNN, J.P., GHENDRIH, P., MONIER-GARBET, P. & POCHEAU, A. 2009 Flow generation and intermittent transport in the scrape-off-layer of the Tore Supra tokamak. *J. Nucl. Mater.* **390**, 368–371.

GARBET, X., PANICO, O., VARENNES, R., GILLOT, C., DIF-PRADALIER, G., SARAZIN, Y., GRANDGIRARD, V., GHENDRIH, P. & VERMARE, L. 2021 Wave trapping and e×b staircases. *Phys. Plasmas* **28**, 042302.

GARBET, X., SARAZIN, Y., BEYER, P., GHENDRIH, P., WALTZ, R.E., OTTAVIANI, M. & BENKADDA, S. 1999 Flux driven turbulence in tokamaks. *Nucl. Fusion* **39**, 2063–2068.

GIANAKON, T.A., KRUGER, S.E. & HEGNA, C.C. 2002 Heuristic closures for numerical simulations of neoclassical tearing modes. *Phys. Plasmas* **9**, 536–547.

GRANDER, F., LOCKER, F.F. & KENDL, A. 2024 Hysteresis in the gyrofluid resistive drift wave turbulence to zonal flow transition. *Phys. Plasmas* **31**, 052301.

GUILLON, P.L. & GÜRCAN, Ö.D. 2025 Phase transition from turbulence to zonal flows in the Hasegawa–Wakatani system. *Phys. Plasmas* **32**, 012306.

GÜRCAN, Ö.D. 2024 Internally driven β-plane plasma turbulence using the Hasegawa–Wakatani system. *arXiv: 2403.09911*.
HAHM, T.S., BEER, M.A., LIN, Z., HAMMETT, G.W., LEE, W.W. & TANG, W.M. 1999 Shearing rate of time-dependent $e \times b$ flow. *Phys. Plasmas* **6**, 922–926.

HALLATSCHEK, K. 2004 Turbulent saturation of tokamak-core zonal flows. *Phys. Rev. Lett.* **93**, 065001.

HASEGAWA, A. & WAKATANI, M. 1983 Plasma edge turbulence. *Phys. Rev. Lett.* **50**, 682–686.

HORNUNG, G., DIF-PRADALIER, G., CLAIRET, F., SARAZIN, Y., SABOT, R., HENNEQUIN, P. & VERDOOLAEGE, G. 2016 $E \times B$ staircases and barrier permeability in magnetised plasmas. *Nucl. Fusion* **57**, 014006.

IDOMURA, Y., URANO, H., AIBA, N. & TOKUDA, S. 2009 Study of ion turbulent transport and profile formations using global gyrokinetic full-f Vlasov simulation. *Nucl. Fusion* **49**, 065029.

IVANOV, P.G., SCHEKOCHIHIN, A.A., DORLAND, W., FIELD, A.R. & PARRA, F.I. 2020 Zonally dominated dynamics and Dimits threshold in curvature-driven ITG turbulence. *J. Plasma Phys.* **86**, 855860502.

KOLMOGOROV, A.N. 1941 The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. *Dokl. Akad. Nauk SSSR* **30**, 301.

KOSUGA, Y., DIAMOND, P.H., DIF-PRADALIER, G. & GÜRCAN, Ö.D. 2014 $E \times B$ shear pattern formation by radial propagation of heat flux waves. *Phys. Plasmas* **21**, 055701.

LABIT, B., FURNO, I., FASOLI, A., DIALLO, A., MÜLLER, S.H., PLYUSHCHEV, G., PODESTÀ, M. & POLI, F.M. 2007 Universal statistical properties of drift-interchange turbulence in TORPEX plasmas. *Phys. Rev. Lett.* **98**, 255002.

LIN, Z., HAHM, T.S., LEE, W.W., TANG, W.M. & WHITE, R.B. 1998 Turbulent transport reduction by zonal flows: massively parallel simulations. *Science* **281**, 1835–1837.

MAHDIZADEH, N., GREINER, F., HAPPEL, T., KENDL, A., RAMISCH, M., SCOTT, B.D. & STROTH, U. 2007 Investigation of the parallel dynamics of drift-wave turbulence in toroidal plasmas. *Plasma Phys. Control. Fusion* **49**, 1005–1017.

MAJDA, A.J., QI, D. & CERFON, A.J. 2018 A flux-balanced fluid model for collisional plasma edge turbulence: model derivation and basic physical features. *Phys. Plasmas* **25**, 102307.

MANZ, P., et al. 2020 Physical mechanism behind and access to the I-mode confinement regime in tokamaks. *Nucl. Fusion* **60**, 096011.

MCDEVITT, C.J., DIAMOND, P.H., GÜRCAN, Ö.D. & HAHM, T.S. 2010 Poloidal rotation and its relation to the potential vorticity flux. *Phys. Plasmas* **17**, 112509.

MCMILLAN, B.F., JOLLIET, S., TRAN, T.M., VILLARD, L., BOTTINO, A. & ANGELINO, P. 2009 Avalanchelike bursts in global gyrokinetic simulations. *Phys. Plasmas* **16**, 022310.

MIKI, K. & DIAMOND, P.H. 2010 Role of the geodesic acoustic mode shearing feedback loop in transport bifurcations and turbulence spreading. *Phys. Plasmas* **17**, 032309.

MOREL, P., GÜRCAN, Ö.D. & BERIONNI, V. 2013 Characterization of predator–prey dynamics, using the evolution of free energy in plasma turbulence. *Plasma Phys. Control. Fusion* **56**, 015002.

NUMATA, R., BALL, R. & DEWAR, R.L. 2007 Bifurcation in electrostatic resistive drift wave turbulence. *Phys. Plasmas* **14**, 102312.

PANICO, O. 2024 Edge turbulence self-organization in fusion plasmas. *PhD thesis*, Institut Polytechnique de Paris, France.

PANICO, O., SARAZIN, Y., HENNEQUIN, P., GÜRCAN, Ö.D., BIGUÉ, R., DIF-PRADALIER, G., GARBET, X., GHENDRIH, P., VARENNES, R. & VERMARE, L. 2025 On the importance of flux-driven turbulence regime to address tokamak plasma edge dynamics. *J. Plasma Phys.* **91**(1), E26.

PARKER, J.B. & KROMMES, J.A. 2013 Zonal flow as pattern formation. *Phys. Plasmas* **20**, 100703.

ROGERS, B.N., DRAKE, J.F. & ZEILER, A. 1998 Phase space of tokamak edge turbulence, the L–H transition, and the formation of the edge pedestal. *Phys. Rev. Lett.* **81**, 4396–4399.

SARAZIN, Y., et al. 2021 Key impact of phase dynamics and diamagnetic drive on Reynolds stress in magnetic fusion plasmas. *Plasma Phys. Control. Fusion* **63**, 064007.

SARAZIN, Y., GARBET, X., GHENDRIH, PH & BENKADDA, S. 2000 Transport due to front propagation in tokamaks. *Phys. Plasmas* **7**, 1085–1088.

SARAZIN, Y. & GHENDRIH, P. 1998 Intermittent particle transport in two-dimensional edge turbulence. *Phys. Plasmas* **5**, 4214–4228.

SCOTT, B. 1997 Three-dimensional computation of drift Alfvén turbulence. *Plasma Phys. Control. Fusion* **39**, 1635–1668.
SCOTT, B.D. 2005 Drift wave versus interchange turbulence in tokamak geometry: linear versus nonlinear mode structure. *Phys. Plasmas* **12**, 062314.

SMOLYAKOV, A.I., DIAMOND, P.H. & MEDVEDEV, M.V. 2000 Role of ion diamagnetic effects in the generation of large scale flows in toroidal ion temperature gradient mode turbulence. *Phys. Plasmas* **7**, 3987–3992.

STAEBLER, G.M., CANDY, J., HOWARD, N.T. & HOLLAND, C. 2016 The role of zonal flows in the saturation of multi-scale gyrokinetic turbulence. *Phys. Plasmas* **23**, 062518.

TERRY, P.W. 2000 Suppression of turbulence and transport by sheared flow. *Rev. Mod. Phys.* **72**, 109–165.

VAN MILLIGEN, B.P., VOLDINER, I., CARRERAS, B.A., GARCÍA, L. & OCHANDO, M.A. 2022 Rational surfaces, flows and radial structure in the TJ-II stellarator. *Nucl. Fusion* **63**, 016027.

VILLARD, L., ANGELINO, P., BOTTINO, A., BRUNNER, S., JOLLIET, S., MCMILLAN, B.F., TRAN, T.M. & VERNAY, T. 2013 Global gyrokinetic ion temperature gradient turbulence simulations of ITER. *Plasma Phys. Control. Fusion* **55**, 074017.

ZHOLOBENKO, W., *et al*. 2024 Tokamak h-mode edge-sol global turbulence simulations with an electromagnetic transcollisional drift-fluid model. arXiv: 2403.10113.

ZHU, H., ZHOU, Y. & DODIN, I.Y. 2020 Theory of the tertiary instability and the dimits shift from reduced drift-wave models. *Phys. Rev. Lett.* **124**, 055002.