2.1 Electron continuity equation

The particle conservation equation, (2.3), evolves the electron density $n$ taking into account the $E\times B$ and diamagnetic drifts, $\boldsymbol {v}_E$ and $\boldsymbol {v}_\star$, the electron parallel current $\boldsymbol {j}_{\parallel,e}$ and the source of particles $S_n$. The polarization drift $\boldsymbol {v}_{{\rm pol}}$, proportional to the mass of the species, is neglected for the electrons

(2.3)\begin{equation} \partial_t n + \boldsymbol{v_E} \boldsymbol{\cdot} \pmb{\nabla} n + n \pmb{\nabla} \boldsymbol{\cdot} \boldsymbol{v_E} + \pmb{\nabla} \boldsymbol{\cdot} (n \boldsymbol{v}_{{\star} e}) - \boldsymbol{\nabla}_\parallel (j_{{\parallel} e}/e) = S_n. \end{equation}

Owing to the large inertia of the ions, the total parallel current is approximated by the electron parallel current: $j_{\parallel e} = - e n v_{\parallel e} \approx j_{\parallel }$. The compressibility terms, $\pmb {\nabla } \boldsymbol {\cdot } \boldsymbol {v}_E$ and $\pmb {\nabla } \boldsymbol {\cdot } (n \boldsymbol {v_{\star e}})$ are kept as derived from (2.1). They are crucial in the linear behaviour of the instabilities as they stabilize the interchange instability at large scale, see § 3. Also, they will prove essential in the energy budget equation § 2.6 and in the nonlinear analysis. In this framework, (2.3) reads as follows:

(2.4)\begin{equation} \partial_t n + \frac{1}{B} \{ \varPhi, n \} - \frac{2 n T_e}{e R B} \left[ \partial_y \left( \frac{e \varPhi}{T_e} \right) - \partial_y \ln \frac{n}{n_0} \right] - \frac{1}{e} \boldsymbol{\nabla}_\parallel j_\parallel = S_n . \end{equation}

With $n_0$ a constant reference density used for normalization and $\varPhi$ the electric potential. The advection of density due to the electric drift is contained in Poisson brackets $\{\varPhi,n\} = (\boldsymbol {\nabla } \varPhi \times \boldsymbol {\nabla } n) \boldsymbol {\cdot } \boldsymbol {e}_\parallel$.

2.2 Charge conservation equation

With the quasi-neutrality assumption, charge conservation reduces to $\pmb {\nabla } \boldsymbol {\cdot } \boldsymbol {j} = 0$. Taking both the diamagnetic current $\boldsymbol {j}_\star = (en \boldsymbol {v}_{\star i} - en \boldsymbol {v}_{\star e})$ and the polarization current $\boldsymbol {j}_{{\rm pol}} = en \boldsymbol {v}_{\rm pol}$ into account, it takes the form

(2.5)\begin{equation} \pmb{\nabla}_\perp \boldsymbol{\cdot} \boldsymbol{j}_\star{+} \pmb{\nabla}_\perp \boldsymbol{\cdot} \boldsymbol{j}_{\rm pol} + \boldsymbol{\nabla}_\parallel j_\parallel = 0 . \end{equation}

Divergence of diamagnetic current can be written directly as

(2.6)\begin{equation} \pmb{\nabla}_\perp \boldsymbol{\cdot} \boldsymbol{j}_\star = \pmb{\nabla}_\perp \boldsymbol{\cdot} ({-}e n \boldsymbol{v}_{{\star} e} + e n \boldsymbol{v}_{{\star} i}) ={-} (1 + \tau) \frac{ 2 T_e}{R B} \partial_y n , \end{equation}

with $\tau = T_i/T_e$ the temperature ratio. The polarization current can be computed from the ion polarization drift, which can be written as

(2.7)\begin{equation} \boldsymbol{v}_{\rm pol} ={-} \frac{m_i}{e B^2} \left[ \partial_t + \boldsymbol{v}_E \boldsymbol{\cdot}\pmb{\nabla} \right] (\boldsymbol{v}_E + \boldsymbol{v_{{\star} i}} ) \times \boldsymbol{B} , \end{equation}

where the parallel advection has been neglected assuming that the parallel gradient length remains small as compared with the transverse ones. This expression of the drift results from the partial cancellation of the advection by the diamagnetic velocity with the anisotropic part of the pressure tensor, which is usually referred to as the diamagnetic cancellation (Hinton & Horton Reference Hinton and Horton1971; Smolyakov Reference Smolyakov1998). Finally, the Boussinesq approximation is used for the polarization current, allowing us to commute the density with the divergence and avoid dealing with time derivatives of different quantities. To keep track of the origin of this density term, we name it $n_\varOmega$ hereafter

(2.8)\begin{align} \pmb{\nabla}_\perp \boldsymbol{\cdot} \boldsymbol{j}_{\rm pol} & \approx{-} \frac{n_\varOmega m_i}{B^2} \partial_t \nabla_\perp^2 \left( \varPhi + \frac{\tau T_e}{e} \ln \frac{n}{n_0} \right) \nonumber\\ & \quad - \frac{n_\varOmega m_i}{B^2} \left( \boldsymbol{v_E} \boldsymbol{\cdot} \pmb{\nabla} \nabla_\perp^2 \varPhi + \pmb{\nabla} \boldsymbol{\cdot} \left( \boldsymbol{v_E} \boldsymbol{\cdot} \pmb{\nabla} \left( \frac{\tau T_e}{e} \boldsymbol{\nabla}_\perp \ln \frac{n}{n_0} \right) \right) \right) . \end{align}

The Boussinesq assumption is considered valid whenever the density gradient length is large with respect to that of the electric potential and is routinely used in fluid models to simplify computation. Note that, depending on the model, the choice is made to commute either a constant density $n_\varOmega =n_0$, such as in Tamain et al. (Reference Tamain, Bufferand, Ciraolo, Colin, Galassi, Ghendrih, Schwander and Serre2016), Dudson & Leddy (Reference Dudson and Leddy2017) and Stegmeir et al. (Reference Stegmeir, Coster, Ross, Maj, Lackner and Poli2018), or to commute the full density $n_\varOmega =n$ (Yu, Krasheninnikov & Guzdar Reference Yu, Krasheninnikov and Guzdar2006; Angus & Umansky Reference Angus and Umansky2014; Ghendrih et al. Reference Ghendrih2022). Although it is usually considered correct in the core, its validity is debated for simulations including the scrape-off layer as the density fluctuations tend to be larger relative to the background density; see dedicated contributions (Yu et al. Reference Yu, Krasheninnikov and Guzdar2006; Angus & Umansky Reference Angus and Umansky2014; Ross, Stegmeir & Coster Reference Ross, Stegmeir and Coster2018). The validity of the Boussinesq assumption can be verified at the end of the simulation by comparing $\boldsymbol {\nabla }_\perp n \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp (\varPhi + \tau \ln (n/n_0))$ and $n \nabla _\perp ^2 (\varPhi + \tau \ln (n/n_0))$. In the simulations performed, the mismatch proved to be of a few per cent. Also, the magnetic field is assumed to commute with the $\boldsymbol {\nabla }$ operator. It is acceptable as the $B$ field is large and decays as $1/R$, which is on much larger scales than the density and electric potential inhomogeneities. Equation (2.5) then reads as follows:

(2.9)\begin{align} & \frac{- n_\varOmega m_i T_e}{e B^2} \partial_t \varOmega - \frac{n_\varOmega m_i}{B^3} \frac{T_e^2}{e^2} \boldsymbol{\nabla}_{{\perp}, i} \left\{ \frac{e \varPhi}{T_e}, \boldsymbol{\nabla}_{{\perp}, i} \left(\frac{e \varPhi}{T_e} + \tau \ln \frac{n}{n_0} \right) \right\} \nonumber\\ & \quad - \frac{2 T_e}{R B} (1 + \tau) \partial_y n + \boldsymbol{\nabla}_\parallel j_\parallel = 0 , \end{align}

where $\boldsymbol {\nabla }_{\perp,i} \{ . , \boldsymbol {\nabla }_{\perp,i} (.) \}$ represents Poisson brackets that can be developed assuming Einstein's notation, $\partial _x \{ ., \partial _x (.) \} + \partial _y \{ ., \partial _y (.) \}$. Here, the generalized vorticity is defined as follows:

(2.10)\begin{equation} \varOmega = \nabla_\perp^2 \left( \frac{e \varPhi}{T_e} + \tau \ln \frac{n}{n_0} \right) . \end{equation}

With $\nabla _\perp ^2 = (\partial _x^2 + \partial _y^2)$ the perpendicular Laplacian.

2.3 Ohm's law

Finally, Ohm's law is derived using the electron parallel momentum conservation equation

(2.11)\begin{equation} n m_e [ \partial_t + (\boldsymbol{v}_E + \boldsymbol{v}_{{\star} i}) \boldsymbol{\cdot} \pmb{\nabla}] v_{{\parallel} e} + T_e \boldsymbol{\nabla}_\parallel n = e n \boldsymbol{\nabla}_\parallel \varPhi + \frac{m_e \nu_{\rm ei}}{e} j_\parallel. \end{equation}

Due to the small electron inertia, the first term proportional to $m_e$ drops. In this case, the generalized Ohm's law reduces to a balance between parallel current, parallel Coulomb force and parallel pressure. One is left with an explicit relationship between the parallel current on the one hand, and density and electric potential on the other

(2.12)\begin{equation} j_\parallel = \frac{e n T_e}{m_e \nu_{\rm ei}} \boldsymbol{\nabla}_\parallel \left( \ln \frac{n}{n_0} - \frac{e \varPhi}{T_e} \right). \end{equation}

The parallel current depends on the electron–ion collision frequency that itself depends directly on the density and temperature: $\nu _{\rm ei} \propto n/T_e^{3/2}$. The temperature is taken constant in the present model but the density has to be made explicit so that we can deal with the pre-factor as a constant. Then, we define $\nu _{{\rm ei},0} = \nu _{\rm ei} n_0 /n_\nu$, the electron–ion collision frequency taken at reference density $n_0$ with the density labelled $n_\nu$ to keep track of its origin.

In this case, the divergence of the parallel current reads as follows:

(2.13)\begin{equation} \boldsymbol{\nabla}_\parallel j_\parallel = \frac{ e n_0 T_e}{m_e \nu_{{\rm ei},0}} \frac{n}{n_\nu} \nabla_\parallel^2 \left( \ln \frac{n}{n_0} - \frac{e \varPhi}{T_e} \right). \end{equation}

The parallel current is then introduced into the electron density and charge conservation equations to close the system.

2.4 Three-dimensional model

The 3-D system of equations then involves the logarithm of the electron density $N=\ln {n / n_0}$ and the generalized vorticity $\varOmega = \nabla _\perp ^2 (\phi + \tau N)$ with the electric potential normalized to the electron temperature $\phi = e\varPhi /T_e$. The dimensionless three-dimensional system can be written as follows:

(2.14)\begin{gather} \partial_t N + \{ \phi, N \} = g \partial_y ( \phi - N)+ \sigma_0 \nabla_\parallel^2 (N - \phi) + D \nabla_\perp^2 N + S_N , \end{gather}

(2.15)\begin{gather}\partial_t \varOmega + \boldsymbol{\nabla}_{{\perp},i} \left\{ \phi, \boldsymbol{\nabla}_{{\perp}, i} ( \phi + \tau N) \right\} ={-}(1 + \tau)g \partial_y N+ \sigma_0 \nabla_\parallel^2 (N - \phi) + \nu \nabla_\perp^2 \varOmega. \end{gather}

Regarding normalizations, time is normalized to the ion cyclotron frequency $\omega _{\rm cs} = (e B)/m_i$ and the lengths to the sound Larmor radius $\rho _s = (m_i c_s)/(eB)$ with $c_s = \sqrt {T_e/m_i}$. The magnetic curvature parameter is defined as $g = (2\rho _s)/R$, with $R$ the major radius of the tokamak. The parallel conductivity is considered constant and is defined as the electron cyclotron frequency to the electron–ion collision frequency taken at $n_0$, $\sigma _0 = \omega _{\rm ce}/\nu _{{\rm ei},0}$. The system is flux driven with a source of particles $S_N$ and the damping of small scales is ensured by the diffusive terms $D$ and $\nu$.

As a last remark, note that (2.14)–(2.15) derive from (2.3) and (2.9) divided by $n$ and $n_\varOmega$, respectively. When doing so, the density dependencies of the parallel current terms should read $(n_0/n_\nu ) \sigma _0$ in (2.14) and $(n_0/n_\nu ) (n/n_\varOmega ) \sigma _0$ in (2.15). The simplifications done in (2.14)–(2.15) are then to replace $(n_0/n_\nu )$ and $(n/n_\varOmega )$ by one. This amounts to performing the Boussinesq approximation with the total density $n_\varOmega =n$, and to neglect the density dependence of the collision frequency $n_\nu =n_0$.Footnote ¹ The latter dependency may be revealed to be particularly important when looking for plasma bifurcations as density increases, such as the density limit. Its effect will be investigated in a forthcoming paper. For the time being, the normalized parallel conductivity is taken as constant and equal to $\sigma _0$.

2.5 Semi-spectral formulation: from three dimensions to one dimension

In order to keep track of the nonlinear dynamics while dealing with a more tractable system, each field is decomposed into a flux-surface-averaged and a fluctuating component (2.16). In the spirit of previous similar models (Sarazin et al. Reference Sarazin, Garbet, Ghendrih and Benkadda2000; Benkadda et al. Reference Benkadda, Beyer, Bian, Figarella, Garcia, Garbet, Ghendrih, Sarazin and Diamond2001; Bian et al. Reference Bian, Benkadda, Garcia, Paulsen and Garbet2003), the fluctuating components are Fourier transformed and projected onto a single parallel and poloidal wave vector $(k_\parallel, k_y)$, so that $(\boldsymbol {\nabla }_\parallel, \partial _y) \longrightarrow {\rm i} ( k_\parallel, k_y)$. Here, c.c. stands for the complex conjugate

(2.16)\begin{equation} \begin{pmatrix} N \\ \phi \end{pmatrix} = \begin{pmatrix} N_{\rm eq} \\ \phi_{\rm eq} \end{pmatrix} (x,t) + \begin{pmatrix} N_k \\ \phi_k \end{pmatrix} (x,t) \exp[{\rm i}(k_y y + k_\parallel z)] + \text{c.c.} \end{equation}

As a caveat, note that the ‘eq’ subscript should not be misunderstood. Here, what is meant by ‘equilibrium’ only refers to the flux-surface average of the considered quantity. This zonal component is not at all assumed to represent any steady state of the system. When discussing $N_{\rm eq}$ and $\phi _{\rm eq}$ in the remainder of this article, we shall loosely refer to them as equilibrium quantities, bearing in mind the caveats mentioned above.

The implications of this choice call for further discussion. In particular, retaining a single poloidal wavenumber $k_y$ implies that, in the time evolution of the fluctuating modes $N_k$ and $\phi _k$, the model cannot consider the nonlinear terms that arise due to mode–mode coupling. Indeed, these terms involve other modes $k_y^\prime \neq k_y$ which are by essence outside the model. One of the consequences is that possible energy (and enstrophy) cascade processes cannot be accounted for. Three important remarks can be made at this point:

1. (i) A refinement of the model would consist in adding a nonlinear saturation mechanism to the fluctuations of the form $\partial _t N_k = \dots -D_{\rm NL} \lvert N_k \rvert ^2 N_k$ with $D_{\rm NL}$ some positive coefficient (Sarazin et al. Reference Sarazin, Garbet, Ghendrih and Benkadda2000). It would account for part of the physics contained in the missing nonlinear interactions, namely nonlinear energy transfer into dissipative scales as one of the routes towards turbulence saturation. We have performed a number of simulations with such a term present, and found that the resulting dynamics was qualitatively analogous to simulations with stronger diffusion. So as to keep the number of free parameters to a minimum, we therefore decided not to include such a nonlinearity in the present study.

2. (ii) Even with $D_{\rm NL}=0$, the model still retains important nonlinearities. The main ones are the turbulent flux and the Reynolds stresses that govern respectively the time evolution of the equilibrium density $N_{\rm eq}$ and poloidal flow $V_{\rm eq}$ profiles (see (2.17)–(2.18) and (2.21)–(2.22)). Note also that, since these radial profiles enter the time dynamics of the fluctuating modes, they result in nonlinear couplings between different radial wave vectors $k_x$ (cf. all terms in (2.19)–(2.20) of the form $F(A_{\rm eq}) G(B_k)$, where $F$ and $G$ stand for linear operators and $A,B \in \{N, V, \varOmega \}$). This latter point is not explicit in the equations since they are not written in the Fourier space in $k_x$ but in the configuration space $x$.

3. (iii) In the absence of an ad hoc nonlinear term in the dynamics of the fluctuating fields (i.e. when taking $D_{\rm NL}=0$), turbulence saturation in Tokam1D therefore relies mainly on two mechanisms: (a) nonlinear transfer of turbulent energy to large-scale flows (ZFs) that do not contribute to transport and (b) the relaxation of mean profiles – as a result of turbulent fluxes – leading to a reduction of the turbulence drive (lower gradients, i.e. thermodynamic forces). The last one is often overlooked since it is absent from gradient-driven models. We argue these two mechanisms play a key role close to marginal stability. As a matter of fact, this regime is likely to be relevant in fusion reactor plasmas. Indeed, given the expected large volume of reactor devices, they will likely be weakly driven – low power density – so that the strong efficiency of turbulent transport – as attested by the stiffness of the experimental temperature profiles in tokamak plasmas (see e.g. Ryter et al. Reference Ryter, Leuterer, Pereverzev, Fahrbach, Stober and Suttrop2000) – should maintain the gradients close to marginality. Note that, since all the modes are kept along the $x$ direction, ZFs can further contribute to turbulence saturation through wave coupling: (I) through turbulence vortex shearing, they can lead to the generation of large $k_x$ wave vectors which are linearly stable (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005); also, (II) they can mediate energy transfer to damped eigenmodes, as highlighted e.g. in Hatch et al. (Reference Hatch, Terry, Jenko, Merz and Nevins2011). Lastly, (III) they can also trap waves (Garbet et al. Reference Garbet, Panico, Varennes, Gillot, Dif-Pradalier, Sarazin, Grandgirard, Ghendrih and Vermare2021), further leading to the possible filamentation of distribution functions in the velocity space and the subsequent collisional dissipation.

The final system of equations involves 4 fields. Two real fields for the equilibrium components (2.17), (2.18), and two complex fields for the fluctuating parts (2.19), (2.20). It is solved using a fourth-order Runge–Kutta scheme in time and centred fourth-order finite differences in the radial direction. Note that the dissipation terms are treated to the second order in the present version of the code, equally limiting the order of the scheme. Appendix D details the numerical scheme and the successful verification tests. Different values for the diffusion and viscosity can be chosen for the equilibrium $(D_0, \nu _0)$ and the fluctuations $(D_1, \nu _1)$

(2.17)\begin{gather} \partial_t N_{\rm eq} ={-} \partial_x \varGamma_{\rm turb} + D_0 \partial_x^2 N_{\rm eq} + S_N \end{gather}

(2.18)\begin{gather}\partial_t V_{\rm eq} ={-} \partial_x \varPi_{\rm RS} + \nu_0 \partial_x^2 V_{\rm eq} - \mu V_{\rm eq} \end{gather}

(2.19)\begin{gather}\partial_t N_k = {\rm i} k_y ( \phi_k \partial_x N_{\rm eq} - V_{\rm eq} N_k ) + {\rm i} g k_y (\phi_k -N_k) + C (\phi_k - N_k ) + D_1 \nabla_\perp^2 N_k \end{gather}

(2.20)\begin{align} \partial_t \varOmega_k & ={-} {\rm i} k_y g (1+\tau) N_k -{\rm i}k_y V_{\rm eq} \varOmega_k + {\rm i} k_y \partial_x[\phi_k \partial_x(V_{\rm eq}+\tau\partial_x N_{\rm eq})] \nonumber\\ & \quad -{\rm i}k_y \partial_x V_{\rm eq} \partial_x(\phi_k+\tau N_k)+ C(\phi_k-N_k) + \nu_1 (\partial_x^2-k_y^2)\varOmega_k , \end{align}

with $\varOmega _k = (\partial _x^2-k_y^2) (\phi _k + \tau N_k)$. The adiabaticity parameter is defined as $C = \sigma _0 k_\parallel ^2$. It is equivalent to the $C_1$ parameter initially derived by Wakatani & Hasegawa (Reference Wakatani and Hasegawa1984). Typically, the parallel wavenumber is estimated considering a connection length, $L_q = 2 {\rm \pi}q R$, which gives $k_\parallel = 2 {\rm \pi}/ L_q = 1/(q R)$. Taking a typical edge plasma of the WEST tokamak, major radius $R = 2.5$ m, minor radius $a = 0.5$ m, magnetic field at separatrix $B_{\rm sep} \approx 3\ {\rm T}$ and considering a density of $n_0 = 10^{-19}\ {\rm m}^{-3}$, a temperature of $T_e \approx 70\ {\rm eV}$ and a safety factor $q_{95} = 5$, the main parameters of the model can be estimated. Their typical values are given in table 1.

Table 1.

Main parameters of the model and their typical values for the WEST tokamak.

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 of $g$ and $C$ values. Nevertheless, one has to keep in mind that (2.12) assumes a large electron–ion collision frequency as compared with the electron inertial term. Therefore, the model loses its validity if $T_e$ is too large or $n$ too small. Also, $C$ decreases with increasing density. This density dependence, although not retained in the current version of Tokam1D, is deemed important for the turbulence-flows energy partition, as highlighted in § 4.2. Finally, we can note that $g=2\rho _s/R=2\rho _\star a/R$, with $\rho _\star =\rho _s/a$. It follows that $g$ is expected to increase for smaller aspect ratio machines operating at the same temperature and magnetic field.

The turbulent particle flux results from the advection of density by the radial $E\times B$ drift

(2.21)\begin{equation} \varGamma_{\rm turb} = \langle n v_{\rm Ex} \rangle ={-} 2k_y \,{\rm Im} (N_k \phi_k^*) ={-}2 k_y \lvert \phi_k \rvert \lvert N_k \rvert \sin \Delta \varphi . \end{equation}

With $^*$ denoting the complex conjugate and $\Delta \varphi = \varphi _k^N - \varphi _k^\phi$ the cross-phase between density and electric potential fluctuations defined as the difference between the phases of the density, $N_k = \lvert N_k \rvert {\rm e}^{{\rm i} \varphi _k^N}$ and of the electrostatic potential. This quantity can be estimated through a linear analysis, see § 3.

The equilibrium flow $V_{\rm eq} = \partial _x \phi _{\rm eq}$ is generated by the Reynolds stress, which is the sum of two components: $E\times B$ and diamagnetic (Smolyakov et al. Reference Smolyakov, Diamond and Medvedev2000; Sarazin et al. Reference Sarazin, Dif-Pradalier, Garbet, Ghendrih, Berger, Gillot, Grandgirard, Obrejan, Varennes and Vermare2021)

(2.22)\begin{equation} \varPi_{\rm RS} = \varPi_E + \varPi_\star = \langle \tilde{v}_{\rm Ey} \tilde{v}_{\rm Ex} \rangle + \langle \tilde{v}_{\rm Ey} \tilde{v}_{{\star} {\rm ix}} \rangle ={-} 2k_y \,{\rm Im}\left[ (\phi_k^* + \tau N_k^*)\partial_x \phi_k \right]. \end{equation}

Finally, a friction $\mu$ is added to the equilibrium velocity equation in order to account for neoclassical damping of the poloidal flow. The code allows for the friction to be estimated through neoclassical assumptions such as in Gianakon, Kruger & Hegna (Reference Gianakon, Kruger and Hegna2002), but we kept it constant in this study. Additionally, the code allows for a relaxation of the poloidal flow towards the radial force balance equilibrium, $\partial _t V_{\rm eq} = \dots - \mu (V_{\rm eq} - V_{\rm eq}^{\rm FB})$, as prescribed in Chôné et al. (Reference Chôné, Beyer, Sarazin, Fuhr, Bourdelle and Benkadda2015). If the toroidal rotation is neglected, the force balanced $E\times B$ velocity reads $V_{\rm eq}^{\rm FB} = v_{\theta } - \partial _x p_i / n$, with $p_i$ the ion pressure and $v_{\theta } = K(\nu _\star, \epsilon ) {\boldsymbol {\nabla }_r T_i}/{e B}$ if one assumes that the poloidal velocity is governed by neoclassical theory. Here, $\nu _\star$ accounts for the collisionality and $\epsilon$ for the inverse aspect ratio. In the case of Tokam1D, due to the isothermal assumption, the velocity from radial force balance would reduce to $V_{\rm eq}^{\rm FB} = - \tau \partial _x N_{\rm eq}$. The coupling of the $E\times B$ flow to the density gradient through radial force balance has already been shown to be able to trigger and sustain edge transport barriers in resistive ballooning mode turbulence simulations (Gianakon et al. Reference Gianakon, Kruger and Hegna2002; Chôné et al. Reference Chôné, Beyer, Sarazin, Fuhr, Bourdelle and Benkadda2015; Giacomin & Ricci Reference Giacomin and Ricci2020). Since our present study aims primarily at investigating how turbulence self-organizes in the presence of self-generated flows, we have decided to drop this radial force balance self-consistency in our study. In this framework, the mean radial electric field in the Tokam1D simulations reported here is only set by the balance between the Reynolds force and the flow damping terms controlled by the friction $\mu$ and viscosity $\nu _0$ coefficients. As a consequence, in particular, the present simulations cannot address possible transitions towards improved confinement regimes such as the H-mode, which has long been known to be associated with a strong modification of the mean radial electric field in the vicinity of the separatrix, in response to the steepening of the pressure profile (Groebner, Burrell & Seraydarian Reference Groebner, Burrell and Seraydarian1990). The impact of the radial force balance will be addressed in a future dedicated work.

2.6 Energetics

Here, we consider the energetics of the system of (2.14), (2.15), which provides a consistency check in order to prevent the appearance of spurious instabilities. Also, partitioning the energy into distinct channels gives insight into the energy transfer between turbulence, flows and profiles in the nonlinear regime. As it should, the model obeys an energy conservation principle, which states that energy is conserved in the absence of a source (particle source $S_N$ here) and sinks (due to dissipative and viscous terms). In other words, linear and nonlinear energy transfer terms are conservative.

To formulate the energy balance equation, we multiply (2.14) and (2.15) by $(1+\tau ) N$ and $(\phi + \tau N)$, respectively. After integrating by parts, and summing the two equations, the interchange terms, proportional to the $g$ parameter, are found to vanish. Such a cancellation requires us to keep compressibility terms in the density continuity equation. The drift-wave terms can be reorganized in the form of a parallel current $j_\parallel = \sigma n \boldsymbol {\nabla }_\parallel (N-\phi )$. The conservation of energy then takes the following form:

(2.23)\begin{equation} \frac{{\rm d}\mathcal{E}_{\rm tot}}{{\rm d}t} = P_\mathcal{E} - D_\mathcal{E}. \end{equation}

With $\mathcal {E}_{\rm tot}$ the total energy, $P_\mathcal {E}$ the production term and $D_\mathcal {E}$ the dissipation term

(2.24)\begin{gather} \mathcal{E}_{\rm tot} = \int E_{\rm tot} \, {\rm d}\mathcal{V} = \int \frac{1}{2} \left\{( 1 + \tau) N^2 + \left[ \boldsymbol{\nabla}_\perp ( \phi + \tau N) \right]^2 \right\} {\rm d}\mathcal{V} , \end{gather}

(2.25)\begin{gather}P_\mathcal{E} = (1+\tau) \int N S_N \, {\rm d}\mathcal{V}, \end{gather}

(2.26)\begin{gather}D_\mathcal{E}= \int \frac{j_\parallel^2}{\sigma_0} \, {\rm d}\mathcal{V} + D (1+\tau) \int (\boldsymbol{\nabla}_\perp N)^2 \, {\rm d}\mathcal{V} + \nu \int \left[\nabla_\perp^2 (\phi + \tau N)\right]^2 \, {\rm d}\mathcal{V}. \end{gather}

With $\int {\rm d}\mathcal {V}$ being the integration on the whole volume $(x,y,z)$. It can be verified that the advection terms, linked to Poisson brackets, vanish upon integration. The production term involves the source of particles injected in the system. The dissipation term contains the parallel plasma resistivity, the dissipation and the viscosity. The total energy $E_{\rm tot}$ involves two terms. The second one scales like the kinetic energy associated with both electric and diamagnetic drifts, while the first one involves both electron and ion pressures. Notice its somewhat peculiar structure already found in Scott (Reference Scott1997), neither proportional to $(1+\tau ^2)N^2$ nor to $(1+\tau )^2 N^2$. Note that the final result, ensuring that the sum is strictly negative, directly results from $\sigma$ depending only on radial direction and time. Therefore, if one chooses to include a non-constant electron–ion collision frequency in the model, it is necessary to consider a flux-surface-averaged quantity.

The energy conservation terms can be decomposed in equilibrium and fluctuating components, following the same method used to derive (2.17), (2.18), (2.19) and (2.20). Noticing that linear terms vanish after the volume integration, the pressure energy can be restricted to

(2.27)\begin{equation} (1+\tau)N^2 \to (1+\tau) (N_{\rm eq}^2 + 2 \lvert N_k \rvert^2) . \end{equation}

The factor 2 in front of $\lvert N_k \rvert ^2$ comes from the definition of (2.16). The generalized vorticity energy reads

(2.28)\begin{align} (\boldsymbol{\nabla}_\perp(\phi + \tau N))^2 & \to \left\{ \left( V_{\rm eq} + \tau \partial_x N_{\rm eq} \right)^2 + 2 \lvert \partial_x \phi_k \rvert^2 + 4 \tau \,{\rm Re}(\partial_x \phi_k \partial_x N_k^*) + 2 \lvert \tau \partial_x N_k \rvert^2 \right. \nonumber\\ & \left.\quad + 2k_y^2 \left[ \lvert\phi_k\vert^2 + \lvert\tau N_k\rvert^2 + 2 \tau \,{\rm Re}(\phi_k N_k^*) \right] \right. \nonumber\\ & \left. \quad - 2 k_y\,{\rm Im} \left[\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^* \right]\vphantom{\left( V_{\rm eq} + \tau \partial_x N_{\rm eq} \right)^2 + 2 \lvert \partial_x \phi_k \rvert^2 + 4 \tau \,{\rm Re}(\partial_x \phi_k \partial_x N_k^*) + 2 \lvert \tau \partial_x N_k \rvert^2} \right\} . \end{align}

The total energy can then be organized into different channels: the equilibrium profiles of density (2.29) and velocity (2.30) on the one hand, and in between fluctuating components on the other hand. The latter can be split into electron and ion (terms proportional to $\tau$) components but for the sake of simplicity, only a turbulent energy that takes all the fluctuations into account (2.31) is considered. Finally, a term corresponding to the interaction between the density and the flows is also obtained (2.32). Its role is still unclear; it remains small in all simulations

(2.29) \begin{align} E_{Neq} &= (1+\tau) N_{eq}^2 + \left( \tau \partial_x N_{eq} \right)^2 \end{align}

(2.30) \begin{align} E_{Veq} &= V_{eq}^2 \end{align}

(2.31) \begin{align} E_{turb} &=2 (1+\tau) \lvert N_k \rvert^2 + 2 \lvert \partial_x \phi_k \rvert^2+ 2 \lvert \tau \partial_x N_k \rvert^2 + 4 \tau {\rm Re}(\partial_x \phi_k\partial_x N_k^*)\nonumber\\ &\quad + 2k_y^2 \left[\lvert\phi_k\vert^2 + \lvert\tau N_k\rvert^2 + 2 \tau {\rm Re}(\phi_k N_k^*)\right]\nonumber\\ &\quad - 2 k_y \Im \left[\left(\phi_k \partial_x \phi_k^*\right) + \tau N_k \partial_x \phi_k^* + \tau \phi_k \partial_x N_k^* + \tau^2 N_k \partial_x N_k^* \right]\end{align}

(2.32)\begin{align} E_{{\rm Neq} - {\rm Veq}} &= 2 \tau V_{\rm eq} N_{\rm eq}. \end{align}

The total energy $E_{\rm tot}$, whose volume integral evolves according to (2.24)–(2.26), can then be recast as

(2.33)\begin{equation} E_{\rm tot} = E_{\rm Neq} + E_{\rm Veq} + E_{\rm turb} + E_{{\rm Neq} - {\rm Veq}} . \end{equation}

It allows one to study the transfer of energy between the density, flows and turbulence. The dynamics between the flow and turbulence energies typically follows a predator–prey behaviour. First, the energy of turbulence increases, until ZFs are created and pump out energy from turbulence. This well-established behaviour has been studied in specific models such as in Malkov & Diamond (Reference Malkov and Diamond2001) and in experiments (Schmitz et al. Reference Schmitz, Zeng, Rhodes, Hillesheim, Doyle, Groebner, Peebles, Burrell and Wang2012). It will be shown that the energy ratio between flows and turbulence greatly changes when scanning the $(g,C)$ parameters of the model: some domains in the parameter space are characterized by large flow to turbulence energy ratio.

3.1 Specific cases: interchange only or CDW only

Setting $\tau = C = 0$ eliminates the drift instability, leaving only the interchange. Consider no dissipation $D=\nu =0$ for simplicity. The dispersion relation then reduces to

(3.9)\begin{equation} \bar{\omega}^2 - \bar{\omega} g k_y - \frac{g k_y}{k_\perp^2} \left[ g k_y - \omega_{*e} \right] = 0 . \end{equation}

The instability develops above the density gradient threshold

(3.10)\begin{equation} \omega_{*e}^{\rm crit} = g k_y \left( 1 + \frac{k_\perp^2}{4} \right) . \end{equation}

In particular, the instability disappears for negative density gradients, recovering the asymmetric nature of the interchange instability with respect to the magnetic field inhomogeneity (Liewer Reference Liewer1985). Also, the largest scales (small $k_\perp$) are found to be destabilized first. Interestingly, the instability threshold increases with $g$. This effect results from the compressibility terms, as will be discussed in the following.

Then, let us consider the purely drift instability case, with finite $C$ and $\tau$ but with $g=0$, assuming the same simplifications. The dispersion relation then reads as follows:

(3.11)\begin{equation} \bar{\omega}^2 + \bar{\omega} \left[ {\rm i}C(1+\tau+k_\perp^{{-}2}) + \tau \omega_{*e} \right] - {\rm i} \frac{\omega_{*e} C}{k_\perp^2} = 0 . \end{equation}

In the cold ion limit, $\tau = 0$, one is left with the Hasegawa–Wakatani system (Wakatani & Hasegawa Reference Wakatani and Hasegawa1984), for which two limits can be distinguished. The case $C \to 0$ corresponds to the hydrodynamic, highly collisional, regime while $C \to + \infty$ corresponds to the adiabatic regime. Keeping only the leading-order terms, the solution in the hydrodynamic case reads

(3.12)\begin{equation} \bar{\omega}_{{\pm}} \approx{\pm} \frac{1+{\rm i}}{\sqrt{2}} \left( \frac{C\omega_{*e}}{k_\perp^2}\right)^{1/2}. \end{equation}

One effectively recovers the resistive drift instability, mostly unstable at large scales ($k_\perp \approx 0$), with the instability growth rate increasing with both the density gradient and the adiabaticity parameter. It exhibits no threshold in the absence of dissipation coefficients, and the instability develops whatever the sign of the density gradient. In the collisionless regime $C\to +\infty$, the system reduces to the Hasegawa–Mima equation and one is left with stable drift waves

(3.13)\begin{equation} \bar{\omega}_+ = \frac{ \omega_{*e}}{1+k_\perp^{2}}, \end{equation}

where $\bar {\omega }_+$ corresponds to the Doppler shifted drift-wave frequency as given in Hasegawa & Mima (Reference Hasegawa and Mima1978). The system is stable and oscillates at the drift frequency. Also notice that, as expected, the turbulent flux (2.21) is strictly zero in this case since the phase shift $\Delta \varphi$ between density and electric potential fluctuations vanishes.

It appears from the asymptotic analysis that the drift-wave instability is stable at both small and large $C$. The interchange instability exhibits a threshold in density gradient and is then stabilized at very large $g$, as further detailed in the following sections.

Also, the two instabilities ($C=0, g=0$) can be distinguished from their linear cross-phase. On the one hand, interchange displays a much larger sine of the cross-phase, typically close to $\Delta \varphi _k \approx {\rm \pi}/2$. This is verified taking $C=0$ and neglecting the compressibility terms in both (3.9) and (3.7). One is left with $F(k, \bar {\omega })=\bar {\omega }/\omega _{*e}$ with $\bar {\omega } = {\rm i} \left ( gk_y\omega _{*e}/k_\perp ^2 \right )^{1/2}$. Then, the phase relation is purely imaginary and the corresponding cross-phase reads $\sin \Delta \varphi _k = -1$. On the other hand, the drift-wave instability exhibits a lower cross-phase which further decreases as $C$ increases. At large $C$, the adiabatic regime forces $N_k \sim \phi _k$. It leads to a vanishing phase shift $\sin \Delta \varphi = 0$ implying a vanishing quasi-linear transport. Taking the limit at large $C$ one can replace $\bar {\omega }$ with the drift-wave solution in (3.7). This leads to ${\rm Im}\,F(k, \bar {\omega }) = 0$ and thus $\sin \Delta \varphi _k = 0$. At low $C$ using (3.12), the phase relation at leading order reads $F(k,\bar {\omega }) \approx - (1+{\rm i}) (C/2k_\perp ^2\omega _{*e})^{1/2}$, which implies $\sin \Delta \varphi _k = -1/\sqrt {2}$.

3.2 General case: coupling CDW and interchange instabilities

Both CDW and interchange instabilities are expected to coexist in the edge of tokamak plasmas. Understanding their coupling is therefore essential. While drift waves are stabilized at both small and large $C$, we show that the interchange instability does not take over in these regimes in the Tokam1D model. Indeed, interchange is also stabilized as the adiabaticity parameter is increased. Recall that the sine of the cross-phase gives an indication of the efficiency of the transport. Thus, predicting its linear value together with the growth rate helps us understand which regimes are expected to yield stronger transport. Finally, we analyse the role of the compressibility terms in the general case, and show that they stabilize the largest scale at very large $g$. The role of $\tau$ is considered in Appendix A. We show that it can be stabilizing or destabilizing depending on the instability at play. Also, it is shown that having a single $k_y$ is detrimental to the study of $\tau$. Therefore, the choice is made to keep $\tau =1$ for the rest of the study.

In the following, the dispersion relation is solved for different cases with fixed equilibrium parameters. Whenever used, the density gradient is fixed to $-N_{\rm eq}'= \rho _s/L_N = 1/100$, $L_N$ being the gradient length. The diffusion parameters are set to $D_1=\nu _1=10^{-2}$. The size of the radial domain is $L_x = 400$ (this constrains $\min (k_x)$). The radial wavenumber is defined as $k_{x} = m ({2 {\rm \pi}}/{L_x})$, with $m$ the radial mode number and $L_x$ the size of the simulation domain.

The growth rate is plotted as a function of $k_x$ and $k_y$ for $C=10^{-3}$ and $g=2 \times 10^{-3}$ in figure 1(a). The white contour corresponds to the linear threshold $\gamma (k_x,k_y)=0$ with these parameters. The white crosses note the positions of the maxima for positive and negative poloidal wavenumbers. In figure 1(b), the effect of $C$ and $g$ on the growth rate is explored. Here, the wavenumbers are fixed: $(k_x, k_y) = (0.06, 0.3)$. Three cases are shown: a CDW only case with $g=0$, an interchange only with $C=0$ and a case ‘both’ with $g=C$.

Figure 1.

(a) Growth rate without equilibrium flows as a function of $k_x$ and $k_y$ for $g=2 \times 10^{-3}$, $C=10^{-3}$. White crosses correspond to maxima in positive and negative $k_y$. White contour denotes the threshold $\gamma (k_x,k_y)=0$. (b) Growth rate for the case collisional drift waves (CDW) $g=0$, interchange (inter) $C=0$ and coupled CDW–inter $C=g$, as a function of $C$ and $g$ for $k_x=2{\rm \pi} /400$ and $k_y=0.3$. Both figures are computed considering: $1/L_N=1/100$, $D_1=\nu _1=10^{-2}$.

In figure 1(a), the growth rate is maximum for the largest radial scale achievable for the system: $k_x = 2{\rm \pi} /L_x$, and finite poloidal scale: $k_y \sim 0.35$. Note that in practice the solution $k_x=0$, although given by the linear analysis, would corresponds to a constant fluctuation along $x$ direction which is imposed to zero by the Dirichlet boundary condition on both boundaries. The growth rate displays a symmetric pattern for positive and negative poloidal and radial wavenumbers. Also, it decays like $-k_\perp ^2$ due to small-scale dissipation governed by the $D$ and $\nu$ coefficients. In figure 1(b), the growth rate of the $(k_x=2{\rm \pi} /400, k_y=0.3)$ mode is plotted as a function of the $g$ and $C$ parameters. The case interchange only leads to a larger growth rate, with $\gamma$ increasing rapidly with $g$. It is abruptly stabilized at $g=10^{-2}$ due to compressibility terms. The CDW are found to be stable, as expected, in both small and large $C$ limits. For the chosen equilibrium parameters, the CDW only case ($g=0$) exhibits a smaller growth rate as compared with interchange. When both instabilities are taken into account, the growth rate stays relatively small. The CDW then proves to have a stabilizing effect on the interchange instability when $C=g$ with the growth rate being barely above the pure drift-wave case.

A more thorough analysis is provided by studying the full parameter space $(C,g)$ while scanning the values of $k_x$ and $k_y$. The radial and poloidal wavenumbers of interest here are those that maximize the growth rate. The result is shown in figure 2 as a function of $C$ and $g$. The white contour corresponds to the threshold $\gamma (g,C)=0$. This figure displays the growth rate $\gamma$, the absolute value of the sine of the cross-phase $\lvert \sin \Delta \varphi \rvert$ and $k_y$ corresponding to the maximum growth rate. The radial wavenumber that maximizes $\gamma$ is always equal to $k_x = 2{\rm \pi} /L_x$.

Figure 2.

Linear analysis of the system without equilibrium flows as a function of $C$ and $g$ for a fixed density gradient. (a) Growth rate $\gamma$. (b) Sine of the cross-phase between density and electric potential fluctuations. (c) Poloidal wavenumber $k_y$ corresponding to the maximum growth rate. The radial and poloidal wavenumbers $(k_x, k_y)$ are chosen such that the growth rate is maximal.

The growth rate governed by the two coupled instabilities, in figure 2(a), is maximal at large $g$ and small $C$: the interchange instability alone leads to a larger growth rate. The drift $C$ parameter acts as a stabilization to the interchange instability. At large $g$, the growth rate is always negative whatever the value of $C$. This is due to the compressibility terms and will be detailed in the next section.

In figure 2(b), the sine of the cross-phase is maximal for interchange dominated cases, at large $g$ for a fixed $C$, consistently with the analytical developments performed in § 3.1. Indeed, the case interchange only is expected to yield $\sin \Delta \varphi = -1$ when the compressibility terms are neglected. In the case of figure 2(b), the compressibility terms do not seem to change this behaviour. The cross-phase decreases with $C$ at fixed $g$. That is also expected from the asymptotic analysis performed earlier. The growth rate is negative at high $C$ whatever the value of $g$. In other word, drift waves dominate at large $C$ while interchange takes over at low $C$. Finally, the most unstable poloidal wavenumber varies from $0.15$ to $0.4$ in the studied parameter domain. The variation remains marginal considering $g$ and $C$ are changed over two decades each.

Remember that Tokam1D considers a constant and unique poloidal wavenumber $k_y$ ($k_y = 0.3$ in the simulations). In turn, this can lead to spurious effects such as the stabilization of the interchange instability at very large gradients. Indeed, since the instability shifts towards higher values of $k_y$ as the gradient increases, the retained $k_y$ may become linearly stable while larger $k_y$ modes get excited. In such regimes, Tokam1D will then consider the system as linearly stable while it is actually still unstable, but at smaller scale (larger $k_y$). However, note that, for physics-relevant parameters, the stabilization does not occur as the density profile stays close to marginality. Moreover, the density starts from a flat profile and slowly builds up with the source. Therefore, there is little chance that the profile suddenly stiffens to reach stabilization. Also, we note that taking a single $k_y$ leads to a stabilization of the interchange instability at very low values of $C$ and $g$, $C<10^{-7}$ and $g<10^{-5}$, which is out of the scanned parameter space. Other than these specific cases, we expect the role of $k_y$ to be marginal on the linear results, and no strong qualitative changes have been observed. A possible upgrade of the reduced model Tokam1D would consist in considering a non-constant poloidal wave vector $k_y$, that would depend on the equilibrium gradients. Such a refinement has not been retained in the present study.

3.3 Compressibility terms stabilize interchange instability at large scales

The compressibility terms, due to the divergence of the electric velocity and diamagnetic flux in the particle balance equation (2.3), involve the magnetic curvature $g \partial _y ( \phi - N)$ (2.14). Those terms stabilize the interchange instability at large scales. Taking the same case as in figure 2 and turning off the compressibility terms, one obtains the growth rate of figure 3(a). In that case the interchange instability is not stabilized and the growth rate diverges with $g$. Additionally, one can look for the minimum gradient needed to destabilize the system figure 3(b). Similarly, due to compressibility terms, this critical gradient increases at large magnetic curvature. Note that there is a threshold for both cases at low $g$ due to the dissipation terms $D$ and $\nu$.

Figure 3.

(a) Growth rate computed using parameters of figure 2 without the compressibility terms. (b) First density gradient $R/L_n$ to destabilize the system as a function of $g$. The red and blue shaded regions correspond to the stable parts of the red and blue curves, respectively.

The stabilization due to compressibility terms can be obtained analytically. A simplified system with, $C=0$, $\tau =0$ and no dissipation or viscosity $D=\nu =0$ leads to the interchange instability. In order to distinguish the compressibility terms coming from the diamagnetic and $E\times B$ terms, the interchange parameter is written as $g_\star$ and $g_E$, respectively, even though both are actually equal to $g$ as they come from the same curvature term. The threshold can be written as

(3.14)\begin{equation} \omega_{*e}^{\rm crit} = g_E k_y \left[ 1 + \frac{g_\star^2k_\perp^2}{4gg_E} \right] , \end{equation}

where it seems that the interchange instability threshold is linearly proportional to compressibility. In the absence of compressibility terms and without dissipation or viscosity, there is no threshold and $\omega _{*e}^{\rm crit}=0$.

4.1 Rich dynamics and flow patterns

In the present model, ZFs are generated nonlinearly by the Reynolds stresses. The action of these sheared flows on the fluctuations then regulates turbulence and therefore the transport. In fact, it is not the amplitude of these flows but whether or not they organize into well-defined sheared layers, leading to corrugations of density or pressure profiles, called staircases, that is the determining factor. It may critically impact the resulting turbulence saturation and therefore the overall quality of confinement. A strong shear may be expected to tilt and elongate turbulent structures, leading to their decorrelation, provided that the shear persists longer than the lifetime of the turbulent structure.

In figure 5 are presented four different velocity patterns as a function of time and radial coordinate $x$ for $C=8 \times 10^{-4}$, $10^{-2}$, $5 \times 10^{-2}$ and $5 \times 10^{-1}$. Each case is taken when the profiles have reached steady state and is plotted for $1.2\times 10^{4}$ time steps, corresponding to $1.2\times 10^{5}$ cyclotron periods $\omega _{\rm cs}^{-1}$. Note that the scale of the colour bars is different for each case.

Figure 5.

Four examples of equilibrium velocity $V_{\rm eq} = - \langle E_r \rangle$ at equilibrium $(\tau > \tau _p)$. Panels show (a) $C=8 \times 10^{-4}$, (b) $C=10^{-2}$, (c) $C=5 \times 10^{-2}$, (d) $C=5 \times 10^{-1}$.

In all simulations, ZFs are active and prove crucial to mitigating turbulence. When the flows are artificially switched off by removing the Reynolds stress drive, a large, system-size radial mode develops. In those cases, the system enters a quasi-periodic regime where a large density gradient builds up and relaxes through a strong transport event. An example of this regime is given in Appendix B. When flows are present, they can be radially structured and stable, as in figure 5(d) or intermittent as in figure 5(c). For the other cases (a,b) the system exhibits intermediate structuration with meandering, splitting and merging events. An example of splitting appears in case (a) at $x=100$ with a fork-like pattern of the flow. Merging occurs in case (b) at $x\approx 220$. The radial structure of the flow also depends on damping. An example is shown in Appendix C, for the case of figure 5(b), where greater flow time stability is obtained when the dissipation coefficients on fluctuations $D_1$ and $\nu _1$ are increased.

The dynamics depends on the adiabaticity parameter: at low $C$, the system develops smooth structures that evolve slowly. When $C$ increases, the dynamics of turbulence and flows is faster. Moreover, diagonal stripes are visible for cases (b,c), associated with ballistic transport events of particles called avalanches. It is known that staircases can act as micro-barriers for these avalanches, efficiently limiting their radial extension (Kosuga et al. Reference Kosuga, Diamond, Dif-Pradalier and Gürcan2014). In general, the flows can either stop the avalanches or, if the avalanches are strong enough, they can travel through the shear layer, perturbing it in the process. In the observed cases, when the flows are radially structured, they always manage to recover their radial structure after an avalanche. Importantly, it is found that the radial size of the ZF is an emergent property independent of the box size $L_x$. For simulations performed at very low $L_x \approx 10$ , a single ZF can fill up all of the simulation domain.

In the confined region of tokamak plasmas, one expects an additional radial electric field that exists alongside this turbulence-driven flow (remember that $V_{\rm eq} \sim - \langle E_r \rangle$ in dimensionless variables). Indeed, the radial force balance relates the mean $E_r$ to the ion pressure gradient and the Lorentz force. This physics is ignored in the present version of Tokam1D, the focus being on the self-regulation of turbulence by self-generated ZFs. Therefore, the flows here are mostly sinusoidal and symmetric around zero.

The diagonal stripes are also visible in the amplitude of the density fluctuations. Figure 6(a) shows the amplitude of the density fluctuations for $C=10^{-2}$, corresponding to the flows in figure 5(b). Density fluctuations are propagating both in and outward at similar speed. The associated turbulent particle flux is always positive, resulting in an outward transport. For the case shown, the typical avalanche length is comparable to the radial size of the velocity structures. The right-hand side features the form of the turbulence in the $(x,y)$ plane for $t=6.03\times 10^6$. It is computed using the complex density fluctuation field and the definition (2.16). On top of the density fluctuation is shown the equilibrium velocity after being coarse grained on a few time and radial points. The tilt of the turbulent structures in the $(x,y)$ plane corresponds to the velocity direction, indicating that flows are efficiently elongating turbulent structures in the poloidal direction. Noticeably, the phase of the fluctuations can exhibit quasi-profile discontinuities, as visible e.g. at $x\sim 160$.

Figure 6.

Density fluctuation for case $C=10^{-2}$. (a) Amplitude of the Fourier component $N_k$ as a function of time and space, the white line notes the snapshot for which is plotted the fluctuation field. (b) Two-dimensional fluctuating density field $\tilde {N}(x,y,t) = N_k(x, t) \exp ({\rm i}k_yy) + \text {c.c.}$ taken at $z=0$ (cf. (2.16)) for $t=6.03\times 10^6$. The black line shows the equilibrium velocity around the same time snapshot.

4.2 Transition from turbulence to flow dominated regime at large $C$

Zonal flows act as a repository for the energy of the system. The more energy is stored inside the ZFs, the less is available for turbulence to generate transport. Numata et al. have shown that there was a collapse of relative energy stored in the flows at low $C$ (Numata et al. Reference Numata, Ball and Dewar2007). Here, we show that, in flux-driven regimes, the density gradient also adapts to the presence of flows and turbulence. Therefore, the transition from a turbulence dominated regime to a flow dominated regime is much smoother. This is in agreement with previous work (Numata et al. Reference Numata, Ball and Dewar2007), provided one accounts for the gradient evolution together with the flows.

Using the energy channels defined (2.29)–(2.32), we can evaluate whether the free energy is captured by turbulence or is stored in the flows. Flows and turbulence energy $E_{\rm Veq}$ and $E_{\rm turb}$ are coarse grained on a few turbulence auto-correlation times. Then, the root-mean-square (r.m.s.) radial profiles are computed

(4.2)\begin{equation} E_{\rm turb}^{\rm rms} (x) = \sqrt{ \langle E_{\rm turb}^2 (x,t) \rangle_t }, \end{equation}

with $\langle \ldots \rangle _t$ a time average. To get a single value per simulation, the radial average of the energy partition is performed, $\langle E_{\rm Veq}^{\rm rms} / (E_{\rm Veq}^{\rm rms} + E_{\rm turb}^{\rm rms}) \rangle _x$. The flow and turbulence energy channels are shown in figure 7(b) as a function of $C$ and their partition ratio is shown in figure 7(a). The colour of the points corresponds to the steady-state density gradient, already shown in figure 4.

Figure 7.

(a) Flow to turbulence energy partition ratio as a function of $C$ with the colour indicating the absolute value of the density gradient in log scale. Error bars represent the standard deviation of the r.m.s. profiles. (b) Corresponding normalized flow and turbulence energies.

The flow to turbulence energy partition ratio increases with the adiabaticity parameter; see figure 7(a). At low $C$ the system is dominated by turbulence with approximately $0.001\,\%$ of the energy stored in the flows. In the adiabatic regime, at large $C$, flows account for $80\,\%$ of the total energy. Consistently, the density gradient increases with $C$. It is illuminating to look at each channel separately. In figure 7(b) the turbulence and velocity channels computed using (4.2) and divided by the total energy are shown. Most of the energy variation is carried by turbulence: it evolves from $3\times 10^{-4}$ of the total energy at low $C$ to $4 \times 10^{-6}$ at large $C$.

On the basis of these observations, we focus on the expected impact of density on the ZF and turbulence dynamics. Although taken constant in the present simulations, the adiabaticity parameter $C=(k_\parallel \rho _s)^2 \omega _{\rm cs}/\nu _{\rm ei}$ should actually scale like $1/N_{\rm eq}$. High density plasmas are then characterized by small $C$ values (assuming constant $k_\parallel \rho _s$). In this case, on the basis of figure 7(a), one expects a low ZF to turbulence energy ratio. This suggests that the turbulent transport should increase when density increases or respectively when $C$ decreases. The trend is qualitatively similar to results in Numata et al. (Reference Numata, Ball and Dewar2007) but not quantitatively. More precisely, the ZF magnitude decreases rather monotonically at large density (collisionality) and does not exhibit the collapse reported in Numata et al. (Reference Numata, Ball and Dewar2007). The difference results from the self-consistent evolution of $N_{\rm eq}$ in flux-driven simulations. In the absence of nonlinear mode–mode coupling , the system has two ways to saturate turbulence: by profile relaxation (transport) and ZF generation. In low $C$ cases, for which the energy stored in ZFs is minimal, the former mechanism is favoured. The richness of the saturation channels permitted by flux-driven simulations thus leads to a less abrupt transition of the system. It should be noted that figure 7(a) is consistent with the ZF collapse reported in Numata et al. (Reference Numata, Ball and Dewar2007) provided one also moves from one density gradient to another when decreasing $C$ (green and red curves in figure 5 of Numata et al. Reference Numata, Ball and Dewar2007). A steep increase of ZF energy has also been found in gyrofluid modified Hasegawa–Wakatani simulations including warm ions (Grander, Locker & Kendl Reference Grander, Locker and Kendl2024). In this latter contribution, mode–mode interactions are retained, leading to a self-consistent turbulent cascade, but the background density gradient is imposed, limiting the liberty of the system to explore the $(C, \partial _x N_{\rm eq})$ parameter space.

A few experiments have reported the strong reduction of ZFs – or proxies for ZFs – when approaching the density limit in L-mode tokamak plasmas (Xu et al. Reference Xu, Carralero, Hidalgo, Jachmich, Manz, Martines, van Milligen, Pedrosa, Ramisch and Shesterikov2011; Hong et al. Reference Hong, Tynan, Diamond, Nie, Guo, Long, Ke, Wu, Yuan and Xu2017; Long et al. Reference Long, Diamond, Xu, Ke, Nie, Li, Wang, Xu and Duan2019). In this context, the ZF collapse at large density in gradient-driven simulations of Hasegawa–Wakatani turbulence (Numata et al. Reference Numata, Ball and Dewar2007) has recently been put forward as a possible explanation for these observations (Hajjar et al. Reference Hajjar, Diamond and Malkov2018; Singh & Diamond Reference Singh and Diamond2022). Our simulations suggest, however, that, in the flux-driven regime relevant to tokamak plasmas, this reduction of ZFs could be much more gradual. Although not discarding the physics as a possible explanation for the issue of density limit in the L-mode, these new results advocate for a renewed exploration of the link between ZF reduction and density limit in a self-consistent flux-driven regime. Whether these findings may explain the reason why ZFs are found to play a secondary role in the density limit mechanism reported in Giacomin et al. (Reference Giacomin, Pau, Ricci, Sauter and Eich2022) remains an interesting yet open question. Last, note that the above analysis has only focused on the impact of density on the adiabaticity parameter, $C$, that governs the transition from the adiabatic to the hydrodynamic regime in Hasegawa-Wakatani turbulence. However, increasing density also increases the friction coefficient acting on the poloidal flow, hence on ZFs ($\mu$ parameter in (2.18)). Gianakon et al. (Reference Gianakon, Kruger and Hegna2002) provide a heuristic expression for this neoclassical coefficient. This density dependency has already been proven to be key to the modelling of Low to High (L-H) confinement bifurcations in flux-driven simulations of resistive ballooning turbulence (Chôné et al. Reference Chôné, Beyer, Sarazin, Fuhr, Bourdelle and Benkadda2014, Reference Chôné, Beyer, Sarazin, Fuhr, Bourdelle and Benkadda2015). This improvement of the model is left for future work.

5.1 Staircase nucleation

Staircases are robust structures. We have found that they always manage to recover after a perturbation. Here, we look at the nucleation of a staircase when restarting a steady-state simulation after having smoothed out the corrugation in $N_{\rm eq}$ and $V_{\rm eq}$ and damped out the fluctuations $N_k$, $\varOmega _k$. The system then has no memory of the past structures. The chosen simulation corresponds to $C=5 \times 10^{-4}$, with the steady-state flows shown figure 5(a). The nucleation of the flow pattern is shown in figure 8(a) together with the dynamics of the density profile in figure 8(b). The black rectangle indicates the time window of the nucleation process which is studied in the following.

Figure 8.

(a) Equilibrium velocity after being smoothed out at $t=3 \times 10^4$. Nucleation is indicated with dotted lines. (b) Evolution of the equilibrium density as compared with its value at the restart: $N_{\rm eq} - N_{\rm eq}(t=3 \times 10^4)$.

Upon restart, at $t=3\times 10^4$ in figure 8, the equilibrium density gradient is above the linear threshold. Therefore, instability and flows appear directly in the whole domain after a short growth time. Qualitatively, a similar staircase structure is recovered. However, there are approximately $7$ structures in figure 8, as compared with the initial $9$ structures.

Focusing on the nucleation of a particular staircase step located at $x=123$, which is indicated with a vertical dashed line, the time evolution of the flows and Reynolds stress at this location are shown in figure 9.

Figure 9.

Energy of the equilibrium flow, total Reynolds stress and its components as a function of time for $x=123$.

An exponential growth of both electric and diamagnetic components of the Reynolds stress is observed. The flows result from a secondary instability (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005) and are governed by the Reynolds stress components. For this simulation, both components of the Reynolds stress are important for the generation of the flows, with the diamagnetic one being slightly larger at initial times.

The dynamics of the phase is essential in the generation of both components. Bearing in mind that fluctuations can be decomposed into amplitude and phase, $\phi _k = \lvert \phi _k \rvert {\rm e}^{{\rm i}\varphi _k^{\phi }}$, the electric component of the Reynolds stress can be recast as

(5.1)\begin{equation} \varPi_E = \langle \tilde{v}_{\rm Ex} \tilde{v}_{\rm Ey} \rangle ={-}2k_y \lvert \phi_k \rvert^2 \partial_x \varphi_k^{\phi}. \end{equation}

The ZF dynamics is governed by the divergence of the Reynolds stress. The separate roles of the phase and of the amplitude clearly appear when computing the logarithmic derivative of $\varPi _E$

(5.2)\begin{equation} \frac{\partial_x \varPi_E}{\varPi_E} ={-}2k_y \frac{2 \partial_x \lvert \phi_k \rvert}{\lvert \phi_k \rvert} + \frac{\partial_x^2 \varphi_k^\phi}{\partial_x \varphi_k^\phi}. \end{equation}

The value of $\varPi _E$ is displayed in figure 10(a) together with its decomposition. The Reynolds force is detailed in figure 10(b).

Figure 10.

Electric Reynolds stress decomposition as a function of time for $x=123$. (a) Reynolds stress. (b) Logarithmic derivative of the Reynolds stress, relating to the Reynolds force. The phase jumps have been removed by ‘unwrapping’ the phase in the post-processing.

In figure 10(a), it appears that both the amplitude and the phase are important for $\varPi _E$ growth. When computing the Reynolds force, see figure 10(b), the phase itself dominates the dynamics. The key role of the phase was already pointed out in Guo & Diamond (Reference Guo and Diamond2016) and in Sarazin et al. (Reference Sarazin, Dif-Pradalier, Garbet, Ghendrih, Berger, Gillot, Grandgirard, Obrejan, Varennes and Vermare2021). It shows that it is essential for a reduced model to retain the complete fluctuations as both amplitude and phase play a role in the generation of flows.

The diamagnetic component of the Reynolds stress, (2.22), can be decomposed in the same way. Considering the conjugate of density fluctuations as $N_k^* = \lvert N_k \rvert {\rm e}^{-{\rm i} \varphi _k^N}$, the tensor reads

(5.3)\begin{align} \varPi_\star & ={-} 2 k_y \tau \,{\rm Im}(N_k^* \partial_x \phi_k) \nonumber\\ & ={-}2 k_y \tau \lvert N_k \rvert \lvert \phi_k \rvert \left[ \partial_x \varphi_k^{\phi} \cos \Delta \varphi + \frac{\partial_x \lvert \phi_k \rvert }{\lvert \phi_k \rvert} \sin \Delta \varphi \right]. \end{align}

With $\Delta \varphi = \varphi _k^{\phi } - \varphi _k^{N}$ the cross-phase between the density and electric potential fluctuations. The first term on the right-hand side is proportional to the electrostatic component of the Reynolds stress while the second relates to the turbulent flux of particles. It writes

(5.4)\begin{equation} \varPi_\star = \tau \frac{\lvert N_k \rvert}{\lvert \phi_k \rvert} \cos \Delta \varphi \varPi_E + \tau \partial_x \left( \log \lvert \phi_k \rvert \right) \varGamma_{\rm turb}. \end{equation}

The correlation between these two components of the total Reynolds stress ultimately depends on the relative weight of the second term with respect to the first. If it is negligible, then the two tensors appear correlated. In this case, the sign of the phase coupling ($\cos \Delta \varphi$) then determines whether $\varPi _\star$ and $\varPi _E$ are in or out of phase. For CDW turbulence, the two tensors appear correlated and in phase ($\mathcal {C}(\varPi _e, \varPi _\star )\approx 0.7$ for $C=3\times 10^{-4}$), which is even stronger at large $C$ ($\mathcal {C}(\varPi _e, \varPi _\star )\approx 0.98$ for $C=8\times 10^{-1}$). This is expected, since the phase between density and electrostatic potential fluctuations decreases approaching the adiabatic limit. The diamagnetic contribution is also larger than the electrostatic contribution with $\varPi _\star > 1.5 \varPi _E$ for $C<10^{-2}$, then slowly decreases towards $\varPi _\star \approx 1.1 \varPi _E$ at large $C$. As a matter of fact, turbulence in the CDW regime is carried mostly by density fluctuations, which results in $\lvert N_k\rvert / \lvert \phi _k \rvert >1$.

5.2 Profile relaxation crucial to self-organization

Profile relaxation is an important saturation mechanism for the turbulence in tokamaks. It is thus crucial to allow this mechanism for studying turbulence saturation and turbulence–ZF interactions in these systems. In order to assess the importance of this mechanism, in this section we artificially prevent any profile relaxation by fixing the density gradient and perform a restart of a steady-state flux-driven (FD) simulation in a gradient-driven (GD) framework. Note that the flux-surface-averaged density gradient does not evolve in time at all in these restarted simulations. This is in marked contrast with the more common GD simulations of the Hasegawa–Wakatani system that one can find in the published literature (see e.g. Majda et al. Reference Majda, Qi and Cerfon2018). Indeed, in these latter cases, the density fluctuations contain a zonal component ($k_y=0$) whose time evolving gradient adds up algebraically to the one (constant in time) set by the background density. Yet, a restoring force prevents this back reaction from balancing the initial gradient – and therefore to annihilate it – so that the self-consistent dynamics of the system gets frustrated. Another difference with respect to more standard simulations – apart from the absence of nonlinear mode–mode coupling in Tokam1D (cf. § 2.5) – is that fluctuations are not assumed to be periodic in $x$ in Tokam1D.

We choose the case $C=5\times 10^{-4}$ and continue the simulation while fixing the gradient. The density profile has been smoothed so that the corrugations are removed, as exemplified in figure 11(b). In figure 11(a), the equilibrium velocity $V_{\rm eq}$ and the density fluctuation amplitude $\lvert N_k \rvert$ are shown as a function of space and time. The first half of the simulation is the steady-state FD simulation. At $t=3\times 10^4$, the equilibrium density profile is smoothed using a Savitzky–Golay filter of order $3$ on a 301 point window size. The second half of the simulation corresponds to the GD regime using the smoothed density profile.

Figure 11.

(a) Equilibrium velocity from GD restart of steady-state simulation at $C=5 \times 10^{-4}$: switch from FD to GD is made at $t=3 \times 10^4$. (b) Density gradient and its radial derivative before and after the GD restart. The density profile before the restart is taken at $t=2.9 \times 10^4$.

Figure 11(a) shows that the flows exhibit roughly the same amplitude in both FD and GD settings. However, their radial structure is very different. The chosen case is radially structured in the steady-state FD case, while no clear layering emerges in the GD case. The flow layers are still present just after the restart but this configuration cannot be sustained in the GD regime. These observations suggest that the system – at least in the regimes explored – must be able to store energy in the density (pressure) profile in addition to the flows in order to develop well-localized flow layers. In other words, the staircase structure requires two critical ingredients: localized shear flow layers in association with large density (pressure) gradients, which appear to reinforce each other. The loss of one of these ingredients prevents the appearance of the whole structure. In the GD regime, the density gradient is frozen in time so that this mechanism is absent. This symbiotic relationship between flow layers and pressure corrugations appears to be instrumental in the staircase dynamics.Footnote ³ The electrostatic potential fluctuations, plotted in figure 12(a), exhibit a slightly larger magnitude in the GD regime as compared with the FD regime. Also, their dynamics looks similar to that of the equilibrium flow in the GD regime, from $t = 3\times 10^4$ onwards. In particular, the radial distance of propagation of avalanches is much larger in the GD regime, sometimes reaching the system size.

Figure 12.

(a) Electric potential fluctuations for the GD restart of the steady-state simulation $C=8$ $10^{-4}$. (b) Radial derivative of electric potential fluctuation phase.

The presence of avalanches in itself may appear surprising. Indeed, avalanches are often understood, by analogy with sandpiles, as resulting from a domino-like effect (Diamond & Hahm Reference Diamond and Hahm1995): a localized strong flux tends to flatten the profile locally and steepens it on both sides due to conservation laws. The strong gradients on both sides then lead to strong local fluxes, further leading to local flattening. The process can repeat over long distances, resulting in the formation of voids and bumps that propagate radially up and down hill, respectively. Local profile relaxation then appears to be key to the whole dynamics. In GD simulations where the density gradient is frozen in time, this mechanism is absent and cannot therefore explain the existence of avalanches on $V_{\rm eq}$. Interestingly, it appears that the electric potential fluctuations magnitude $|\phi _k|$ and phase gradient $\partial _x\varphi _k^\phi$ exhibit an avalanche dynamics, as can be seen in figure 12. They both govern the electric component of the Reynolds stress, (5.1). It is so far unclear whether one or the other – or both – is the main drive for the dynamics of $V_{\rm eq}$, the other quantity being slaved to it via the back reaction of the flow shear. Note that both terms, the phase gradient and the fluctuation magnitude, are contributors to the turbulent energy. Therefore, discriminating one scenario with respect to the other is not obvious solely from an energetics standpoint. This remains an open issue for future works.