Investigation of the interaction between competing types of nondiffusive transport in drift wave turbulence

The sustainable confinement of fusion plasmas in a tokamak concept has been a topic of active research for many years now. Inside the separatrix, radial transport in fusion plasmas typically dictates plasma confinement, stability, and plasma-wall interactions. Radial transport in most fusion plasmas has shown to be dominated by turbulence, resulting in transported quantities that mix faster than molecular diffusion, and a broadband spectrum of fluctuations.¹ Hence, there is a need to understand turbulent transport in tokamak plasmas in order to achieve and maintain fusion plasmas for longer durations in future machines.

In the past few decades, the dynamics of turbulent transport in tokamak plasmas have been shown to be much richer than a mere increase in transport coefficients. In particular, the nature of radial transport itself may often become non-diffusive.² Numerical simulations have shown that, in near-marginal regimes in which the separation of timescales between fluctuations and profile modification is importantly reduced, radial transport can be endowed with superdiffusive features in which radial transport takes place via avalanching.^3–6 On the other hand, numerical simulations have also shown that radial transport through regions with strong, radially sheared zonal flows can become subdiffusive.⁷ The physical ingredients that seem to be responsible for these behaviors are quite different. In the superdiffusive cases, a proper profile evolution in the presence of a local instability threshold is essential to get the correct transport dynamics. In the subdiffusive one, the self-consistent generation and evolution of the zonal flows are needed. However, making numerical simulations in which both profile evolution and zonal flow dynamics are done self-consistently and simultaneously is very difficult because of the huge computational requirements. As a result, most numerical simulations of tokamak turbulence traditionally use fixed background profiles, thus including only zonal flow dynamics. In fact, it was in this type of simulations that the subdiffusive nature of radial transport across zonal flows was identified.⁷ Much less frequently, simulations are done that evolve background profiles self-consistently with the turbulence, using what is known as a flux-driven setup. Although flux-driven simulations are becoming more common due to the availability of more powerful supercomputers, even in gyrokinetics,^8,9 still just a handful of these runs have focused on studying the non-diffusive nature of radial transport. Regretfully, zonal flow dynamics have often been played down or even removed to simplify the calculations, as was the case of the simulations in which the superdiffusive nature of radial transport was first identified in plasma turbulence.⁴ 

The importance of understanding the complete interaction between profile modification and zonal flow dynamics and assessing the resulting transport dynamics depending on their relative strength is however self-evident when considering that the ITER tokamak, due to its larger temperatures, will probably operate frequently in near-marginal conditions, where profiles should be expected to vary in timescales not that separate from those of turbulence. In addition, it is also expected that self-consistent zonal flows will lend a hand in keeping turbulent transport under controlled in ITER. Flux-driven, gyro-kinetic simulations in realistic geometries will be needed to undertake this type of studies. However, the past experience has often shown that gyrokinetic simulations are not only very expensive but also very complicated to interpret due to the large number of elements included, both physical (geometry, physics included, etc.) and numerical (algorithms used, simplifications to speed up calculations, etc.). Thus, there is much to be gained by having a first glimpse, in the context of nondiffusive transport, at how the interaction between profile evolution and self-generated flows may look like in a simplified setup, where the relevant ingredients can be easily isolated from other complications and where the knobs that allow tuning its relative importance could be more easily identified.

This paper intends to do precisely this. It carries out a study of this interaction in a simple slab drift-wave turbulence model, extensively studied in the past in a fixed-gradient setup,^10,11 to which simultaneous and self-consistent profile evolution and zonal flow dynamics have been added. As will be shown in what follows, quite meaningful results are readily available from this model. Simply by modifying the parameters that define it, different transport regimes are accessible in which the mutual interaction between flows and profile relaxation changes, and this have allowed us to understand better what determines the final nature of radial transport and also how to modify it. We expect these results to be extremely helpful to shed light on the underlying physics and facilitate the interpretation of future results from much more complex, flux-driven, fluid, and gyrokinetic simulations.

This paper is organized as follows: Sec. II describes both the dissipative trapped electron mode (DTEM) turbulence model that will be used and how profile relaxation and zonal flow dynamics have been added to it; also, the free knobs available to us in order control to the interaction between them will be discussed; Sec. III describes the diagnostics that will be used to characterize the nature of the turbulent transport in each of the different simulations analyzed. In particular, the meaning of the transport exponent H, also known as the Hurst exponent,¹² will be explained. Sections IV and V contain the main results of this study. First, we describe in Sec. IV the main features of the steady-states reached by several simulations for a representative collection of free parameter values; then, in Sec. V, the nature of transport is characterized by determining the transport exponent H of each simulation. A coherent explanation of all these results will be given in Sec. VI, with special focus on pointing out the relevant physics in each case. Finally, this paper concludes by speculating on how our results may translate to more realistic tokamak plasmas in Sec. VII.

In this work, a drift-wave model is used for the turbulence due to its ubiquitous nature, sometimes categorized as “universal instabilities,” as they occur in most confined plasma configurations. As is known, any non-uniform density plasma maintained by a strong magnetic field is susceptible to drift-wave instabilities.¹³ The drift-wave instabilities access the thermal energy as the plasma expands across the magnetic field, and many of their characteristics, including its turbulence spectrum and nonlinear cascades,^11,14 have been quantified numerically in the past using fixed-background setups. The DTEM forms the basic structure of our model. It is derived from the two-fluid plasma equations by assuming cold ions, trapped electrons, and quasineutrality.^10,15–17 The geometry used here is that of a doubly periodic two-dimensional slab, perpendicular to a constant magnetic field, with $x ∈ [ 0 , 1 ]$ as the radial direction and $y ∈ [ 0 , 1 ]$ as the poloidal direction.

The evolution equations of our model are

where $ϕ$ is the fluctuating potential, n is the fluctuating density, and P is the background density. The first two equations are very similar to the standard DTEM model^11,14 but include an additional dependence on the background profile P via the nonlinear function $R n l ( n , ϕ , P )$⁠, which we will discuss later. The third equation represents the profile evolution and includes a nonuniform (in radius) source term S, also discussed later, in order to implement a flux-driven setup. The coupling with the fluctuations happens through the Lagrangian derivative, $d / d t = ∂ t + u · ∇ ⊥$⁠, where $u = C s ρ s z ̂ × ∇ ⊥ ϕ$⁠, the usual turbulent E × B drift. The meaning of the coefficients in the model is, otherwise, pretty standard. $ρ s = ( k B T e / e B ) / C s$ is the ion gyroradius, $C s = k B T e / m i$ is the ion sound speed, ϵ is the inverse machine aspect ratio that defines the trapped electron fraction in velocity space $v ∥ < ϵ v$⁠,^17,μ is the viscosity coefficient, ν is the electron collisional relaxation due to trapping and detrapping, and $ν e f f = ν / ϵ$⁠. The term $ξ = ( 1 + α η e )$ contains the instability criterion, where $η e [ = 2 ]$ is the ratio between the electron temperature gradient and the electron density gradient and $α = 3 / 2$ for the instability criterion for destabilization of DTEM modes by electron collision.^10,13 Coordinates are normalized such that $t → t Ω i , x → x / 10 ρ s$⁠, and $y → y / 10 ρ s$⁠. Finally, the parameter ν specifies the toroidal force balance that translates as a phase lag between n and $ϕ$⁠.

We proceed to discuss first the source term in the background profile evolution equation. The periodic (x and y) source S used (see Fig. 1) has been chosen to introduce the radial inhomogeneity needed to generate zonal flows, as well as to provide the means to establish a flux-driven, evolving background profile. It is defined as

where G(x) is a narrow, normalized Gaussian distribution and S₀ represents an injection rate. Thus, injection is centered around $x source [ = 0.25 ]$ and extraction at $x sink [ = 0.75 ]$⁠. At the steady state, the net amount lost at the sink will balance, on average, what is gained at the source. The background profile will evolve consistently with this balance. Since the source only depends on x, the gradient of the profile is mostly directed along the radial direction, with weaker poloidal dependencies induced by the turbulence. It must also be noted that, since the background profile equation is a convective equation involving energy injection at primarily low-k wavenumbers, a diffusion term proportional to D_P has been added to prevent the formation of extremely steep gradients at high k.

Next, we discuss the two nonlinearities in the DTEM model: the so-called E × B nonlinearity, which results from the advection of the fluctuating density $u · ∇ ⊥ n = C s ρ s z ̂ × ∇ ⊥ ϕ · ∇ ⊥ n$⁠, and the polarization drift nonlinearity, which arises from the advection on the normalized vorticity $u · ∇ ⊥ ( ∇ ⊥ 2 ϕ ) = C s ρ s 3 z ̂ × ∇ ⊥ ϕ · ∇ ⊥ ( ∇ ⊥ 2 ϕ )$⁠. The main characteristics of these two nonlinearities have been extensively studied numerically in a fixed-background context.^10,14,18 The E × B nonlinearity is dominant at small k wave numbers and provides a non-local cascade of energy towards large k. In contrast, the polarization drift nonlinearity is dominant at large k and causes a cascade of energy towards small wave numbers.¹⁰ The nonlinearities do however act to transport momentum in a conservative way, therefore without changing the total momentum.

Finally, we describe the nonlinear function $R n l ( n , ϕ , P )$ that sets up the coupling between the background profile equation and the local turbulence. It is defined as

To understand the meaning of this equation, it is better to focus first on just its first line. The main ingredient is the function $g ( L s , z − 1 )$⁠, which provides the threshold condition on the background profile in order to locally drive the instability, and this is a key element for the establishment of superdiffusive transport dynamics.^3,19 It is a function of the gradient of some arbitrary quantity, $s = s ( z )$⁠, that we express as $L s , z − 1 = ∂ z s / s 0$ [s₀ is just a normalization constant], and this becomes non-zero (which then makes R_nl equal to its fixed-background value for the DTEM model) only above a prescribed threshold value for this gradient. More specifically, g is constructed as a symmetric combination of hyperbolic tangent functions (Fig. 1, right frame)

where $L c , s , z − 1$ is a prescribed critical value for the gradient of field s in the z direction and κ denotes the steepness of the hyperbolic tangent function. [The value of the steepness κ has been fixed at κ = 20 in order to model a step function more closely.] The symmetric combination of tangent functions ensures that instability is independent of the sign of the local gradient, only on its magnitude.

Therefore, the threshold function g appearing in the first line of Eq. (3) introduces a critical radial gradient in the poloidally averaged background profile, $⟨ P ⟩ y$⁠. On the other hand, the factor that multiplies g, $L ⟨ P ⟩ y , x − 1 ∂ y ϕ$⁠, represents the coupling between the local turbulence and the background profile. It is obtained by evaluating the fluctuating turbulent advection of the mean profile,

that, when assuming a very fast equilibration along the field lines (here, along the y direction), reduces to

The second line in Eq. (3) has been added to allow us to explore variations of the threshold term just described. In particular, we are interested in reproducing possible tokamak-relevant situations in which the equilibration along the field lines may not be so fast compared with turbulent timescales; in those cases, magnetic surface equilibration (that is, poloidal and toroidal equilibration) of background profiles should not be expected. We must however introduce this possibility in the model in a rather ad-hoc way since the magnetic field is here perpendicular to the xy plane, and there is no parallel transport. For that reason, it is done by considering a linear combination of two threshold conditions whose relative importance is weighted by the factor $f d ∈ [ 0 , 1 ]$⁠. The first condition, given by the first line in Eq. (3) just described, includes a threshold condition on the radial (i.e., along x) gradient of the poloidally (along y) averaged background profile, $⟨ P ⟩ y = ℓ y − 1 ∫ 0 1 P ( x , y ) d y$⁠. Thus, it is somewhat analogous to a tokamak situation in which fast parallel transport allows for equilibration of any poloidal inhomogeneities. The second line of Eq. (3) introduces the second threshold condition as a function of the local poloidal (in y) and radial (in x) gradients instead. Thus, this term would represent situations in which the lack of sufficiently fast parallel transport makes relevant the response of local turbulence to local gradients. These local gradients modify the local drift velocities that perform the advection on P, thus resulting in a more inhomogeneous, anisotropic drive for the turbulence, leading to larger Reynolds stresses to drive sheared flows. The local threshold condition also adds an additional mechanism to equilibrate gradients via parallel motion, which is absent in the parallel equilibration term proportional to $n − ϕ$ in the potential equation.

In summary, the set of equations of the DTEM model described in Eq. (1) contains all the elements needed to produce both near-marginal transport (through the profile evolution equation and the threshold condition included in the evolution equations for the fluctuations) and transport across zonal flows (that will be generated via the Reynolds stress represented by the two nonlinearities of the DTEM model) simultaneously and self-consistently. Table I summarizes the added parameters as part of the flux-driven setup. The free parameters of the model are the drive strength, S₀, the parallel equilibration factor f_d, and the critical gradient threshold $L c , P , x − 1$⁠. By changing their values, one can increase [or decrease] the intensity of the zonal flows and the distance of the background profile from marginality. In what follows, we will describe several simulations (see, for example, Fig. 2) that explore these parameter variations and investigate the resulting changes in the nature of radial transport. This analysis will allow us to understand better the interplay between self-generated flows and near-marginal turbulent transport.

To conclude this section, we provide some details about the numerical scheme used to solve Eq. (1). The spatial domain considered is a doubly periodic grid of 256 × 256 nodes in the Fourier space spanned by k_x and k_y. The scheme used is a standard spectral one, properly modified to avoid any aliasing problems, that uses the pseudo-spectral method to deal with nonlinearities. The temporal integration is done implicitly, using a scaled preconditioned Generalized Minimal Residual (GMRES) solver that combines well-established integration schemes.²⁰ Parallelization is achieved by using Message Passing Interface (MPI) and by taking advantage of parallel Fast Fourier Transform (FFT) routines²¹ and other parallel numerical integration routines. All simulations have been initialized with random phases for all Fourier harmonics, and they have been advanced in time until a suitable quasi-steady state, with approximate balance between drive and losses, is established. As an illustration, some snapshots of the vorticity field obtained at the steady state of representative simulations used in this paper are shown in Fig. 2. In the figure, the position in the radius of sink and source is shown by means of two vertical, dashed red lines.

There are many ways to define the characteristics of the nature of transport in a system.^22,23 Traditionally, transport of any conserved quantity P is termed diffusive if it can be described by an equation of the type

where χ is a transport coefficient. This equation implies that the local flux of the transported quantity is $Γ P = − χ ∂ x P$⁠, thus pointing down the gradient, which is known as Fick's law. By extension, some authors define nondiffusive transport as any situation in which the evolution of P follows instead a transport equation of the form:

where $β ∈ ( 0 , 1 ]$ and $α ∈ ( 0 , 2 ]$ are known as fractional transport exponents, and $χ α , β$ is the fractional transport coefficient. The operators that appear in this equation are fractional derivatives, and they provide smooth interpolants in between integer derivatives.²⁴ In contrast to the integer derivatives, which are local operators, fractional derivatives are integro-differential operators that are non-local in their variable of definition (either space or time). Thus, the type of transport that can be captured by Eq. (8) may be intrinsically non-local (if $0 < α < 2$ and non-Markovian if $0 < β < 1$⁠). However, the exponent of most interest to us is the exponent $H = β / α$⁠. Transport is termed superdiffusive when H > 1∕2, diffusive if H = 1/2 (even if $β ≠ 1$ and $α ≠ 2$⁠, as in the usual diffusive equation), and subdiffusive when H < 1∕2. The reason for these names is that any population of particles, whose transport is governed by Eq. (8) and that are initially localized in x, will spread faster (if H > 1∕2) or slower (if H < 1∕2) than its diffusive counterpart.

As with the classical diffusion equation (Eq. (7)), which can be derived from “microscopic considerations” that assume Gaussian fluctuations and lack of memory, the fractional transport equation (Eq. (8)) can be obtained by assuming “microscopic” Lévy distributed fluctuations and self-similar memory. It is well known that the classical diffusion equation can be obtained as the long-term, long-distance limit of either a continuous time random walk (CTRW) in which particles are advanced with Gaussian-distributed step-sizes, separated by exponentially distributed waiting times. Also, the diffusion equation can be obtained from the “microscopic” Langevin equation, which expresses the position of each particle being transported as

In the same long-term, long-distance limit, by assuming a random forcing ξ with Gaussian statistics and zero time correlation except at zero lag, the classical diffusion equation is recovered.

Similarly, the fractional transport equation can be derived from a [still badly called random] CTRW²⁵ in which steps are distributed according to a symmetric α-Lévy distribution with tail index $α ∈ ( 0 , 2 )$ and waiting-times distributed according to an extremal β-Lévy distribution with tail index $β ∈ ( 0 , 1 )$⁠.^22,23,26 It can also be derived as the long-term, long-distance limit of the generalized Langevin equation^27,28

that assumes a non-random forcing with symmetric α-Lévy statistics and a correlation in time characterized with a Hurst exponent $H ∈ ( 0 , 1 )$⁠.¹² In that case, $β = α H$ in the transport equation.

The connections of Eq. (8) with “microscopic formulations” can be exploited to come up with methods to determine the fractional exponents α, β, and H in practical situations and thus to provide ways to characterize the nature of transport. Among different methods available, in this paper, we will focus on one that exploits the connection with Eq. (10), permitting us to easily determine the exponent we are interested in, H, by following the trajectories of massless tracer particles as they are advected by the turbulence. That is, by integrating in time their velocity, that is given by

since the advection in our model is done by the turbulent fluctuating E × B velocity.

In order to estimate the H exponent, we benefit from the fact that it represents the self-similarity exponent of the process described by Eq. (10). Therefore, it can also be obtained also as the correlation (or Hurst) exponent of its derivative or time series of the increments of the process or, in discrete form, of its velocity series along the Lagrangian trajectory. The method we have chosen to determine this correlation exponent is the well-known R/S technique, which has been reliably used for more than sixty years.¹² Given a velocity series of length N represented as ${ V k : 1 ≤ k ≤ N }$⁠, one just needs to construct the so-called rescaled range

where $X k = V 1 + ⋯ + V k − k V ¯ n$ for $k ≤ n$⁠. The mean is defined as $V ¯ n = ∑ i = 1 n V i / n$⁠. The s-variance is defined as $σ n = ∑ i = 1 n ( V i − V ¯ n ) s / ( n − 1 )$⁠, with $s < α$⁠, being $α < 2$ the tail-exponent characterizing the velocity statistics (or, if Gaussian-distributed, α = 2). It then happens that, if the signal is correlated with Hurst exponent H (and therefore, its integrated path is self-similar with the same exponent), one finds that

from which the exponent is readily obtained. It is fair to say that the R/S method has been criticized in the literature because it tends to somewhat overestimate exponents (for instance, R/S tends to yield $H ∼ 0.55$ for random signals instead of 0.5), but it is extremely resilient to both noise and periodic perturbations,²⁹ which is why it is our method of choice. It is also worth to note that the statistics of the determination of H are greatly improved by averaging the rescaled range for time n over all non-overlapping segments of size n in which the full time series can be broken, procedure that we have extensively used in this work.

In this section, we will describe the main features of the steady-state background profiles obtained for different sets of simulations examined. Each simulation, once suitable values for the local critical gradient $L c , P , x − 1$ (⁠ $L c , P , y − 1 = 0$ is always used, so a poloidal drive for turbulent transport always exist), the source injection rate, S₀, and the parallel equilibration fraction, f_d, have been prescribed, was advanced in time until an steady-state was reached in which sources and sinks balanced and the background profile P fluctuated around a well-defined average shape. As diagnostics, several (spatially integrated) quantities of interest were monitored, including $W P ∝ ∫ 0 1 ∫ 0 1 | P | 2 d x d y$ and $W turb ∝ ∫ 0 1 ∫ 0 1 ( | n | 2 + | ϕ | 2 ) d x d y$⁠, which provide proxies for the total background and fluctuating energies, respectively. The arrival to the steady-state was revealed by having $⟨ W P ⟩ t ∼ constant$ and $⟨ W turb ⟩ t ∼ constant$⁠, as well as $⟨ ∂ x P ⟩ y , t ∼ constant$ (i.e., a well-defined mean background profile), over many eddy turnover times. The level of turbulent activity was monitored using the quantities $δ W turb / ⟨ W turb ⟩ t$ and $δ W P / ⟨ W P ⟩ t$⁠, where $( δ W turb ) 2 = ⟨ ( W turb − ⟨ W turb ⟩ t ) 2 ⟩ t$ and $( δ W P ) 2 = ⟨ ( W P − ⟨ W P ⟩ t ) 2 ⟩ t$⁠, and with $⟨ · ⟩ t$ standing for time averaging. These quantities are related to the overall amplitude of the turbulence and its effect on distributing the free energy stored in the background profile, respectively. It must be kept in mind that, in systems with both linear and nonlinear terms, relatively large values of $δ W turb / ⟨ W turb ⟩ t$ reveal bursty turbulent activity, while large $δ W P / ⟨ W P ⟩ t$ values are associated with very active background profile modifications, both effects being associated with the action of nonlinear terms. The time-average shape of the gradient of the (poloidally averaged) background profile, $⟨ ∂ x P ⟩ y , t$⁠, is coupled to the values of these quantities as well. Thus, in cases in which $δ W turb / ⟨ W turb ⟩ t$ is large, which implies large levels of turbulence-induced transport, smaller $⟨ ∂ x P ⟩ y , t$ values are obtained for the same source. On the other hand, large $δ W P / ⟨ W P ⟩ t$ values are often indicative of near-marginality, that is, of profiles locally wandering around the local instability threshold values. Finally, small values of both $δ W turb / ⟨ W turb ⟩ t$ and $δ W P / ⟨ W P ⟩ t$ are usually a reflection of a dominant contribution from the linear terms, and larger values of the gradient of the background profile are then expected.

In order to examine the change in transport dynamics in a context of competing zonal flow development and turbulent profile modifications, we have run three different sets of simulations, in which just one of the three available free parameters is varied while keeping the other two fixed. The explored range of parameters is shown in Table II. By increasing from zero the critical gradient threshold, $L c , P , x − 1$⁠, we were able to explore states in which the dominance of self-generated flows is gradually replaced by more prominent radial turbulent relaxations (or avalanches). In addition, scenarios in which the source injection rate S₀ was varied test the resiliency of the avalanche relaxation dynamics. Finally, the cases in which the parallel equilibration f_D was increased from zero assess the importance of the parallel equilibration dynamics on the transport relaxation process.

The effect on the steady-state radial gradient of the background profile, $⟨ ∂ x P ⟩ y , t$⁠, of using different values of $L c , P , x − 1$ for a fixed source $S 0 = 5$ is shown in Fig. 3. In all cases, f_d = 0, which means that parallel equilibration is not being considered. The steady-state profile in the (almost) absence of turbulence is shown in yellow, as obtained by considering the regime where the turbulent advection is small compared to the diffusion term $D P ≫ u ρ s$⁠. Since the diffusive flux is proportional to the local gradient and the whole source at $x source = 0.25$ had to be transported diffusively to the sink location at $x sink = 0.75$⁠, this case exhibits the largest gradients of all in the central region.

The case run with $L c , P , x − 1 = 0$ (no critical gradient, shown in magenta circles) corresponds to the case in which the profile is always unstable and turbulence is always active. This regime is known as one of the supermarginal turbulences. On average, $⟨ ∂ x P ⟩ y , t$ sits at a much lower value than the diffusion dominated profile in between the sink and the source. This situation corresponds, as we mentioned earlier, to significant values for $δ W turb / ⟨ W turb ⟩ t$ in the central region (or the order of $∼ 5 %$⁠), a signal for significant levels of turbulent transport. On the other hand, $δ W P / ⟨ W P ⟩ t$ is much smaller which, as we will see soon, implies a modest level of profile modification (⁠ $δ W P / ⟨ W P ⟩ t ≃ 0.9 %$⁠). This level of fluctuations is sufficient to establish significant radial turbulent fluxes that can transport P from the source to the sink while maintaining smaller gradient values than in the diffusive case, as clearly seen in Fig. 3. The situation is quite different, however, near the source and sink locations. The local inhomogeneity and anisotropy induced by the shape of the source term (Eq. (2)) in those regions drive strong radially sheared, poloidal flows (see Fig. 4) that reduce turbulent fluctuations in those regions to a minimum. Therefore, the radial turbulent transport in these regions is very small, which causes the gradients to increase to diffusive levels in order to maintain the flux of excess free energy through the diffusive channel.

We discuss the next cases with a finite critical gradient threshold for the same source S₀. Three different values have been studied shown in Fig. 3 using blue diamonds (⁠ $L c , P , x − 1 = 0.4$⁠), brown triangles (⁠ $L c , P , x − 1 = 0.8$⁠), and green asterisk (⁠ $L c , P , x − 1 = 1$⁠), respectively. Focusing first on the central region between the source and sink, it seems clear that increasing the value of the threshold initially leaves the average background gradient $⟨ ∂ x P ⟩ y , t$ quite unchanged (⁠ $δ W turb / ⟨ W turb ⟩ t ≃ 7 %$ and $δ W P / ⟨ W P ⟩ t ≃ 0.9 %$ for $L c , P , x − 1 = 0.4$⁠). However, for the largest two threshold values, the background gradient has decreased quite a bit from the $L c , P , x − 1 = 0$ case (indeed, $⟨ W P ( L c , P , x − 1 = 0 ) ⟩ t / ⟨ W P ( L c , P , x − 1 = 1 ) ⟩ t ∼ 25$ and $δ W P / ⟨ W P ⟩ t ≃ 20 %$⁠). This reduction in the total free gradient energy is due to more vigorous turbulent relaxations (indeed, $δ W turb / ⟨ W turb ⟩ t ≃ 10 %$ for $L c , P , x − 1 = 1$⁠.) associated with significant turbulent profile modification, which is the first telltale evidence of transport having reached near-marginality. This is also apparent from Fig. 3, where the (properly normalized) critical threshold value or $L c , P , x − 1 = 1$ has been indicated in red. Clearly, the average gradient shows regions above and below that threshold line, as expected in near-marginal conditions. In contrast, for the smallest thresholds, profiles are above marginal. Therefore, its gradient is set by the source and the level of fluctuations achieved, instead of by the local gradient threshold.

A characteristic property of near-marginal turbulent steady-states is profile resilience or, in other words, the lack of sensitivity of the profile share to the strength of the drive, being instead determined mostly by the local threshold value. To investigate whether our model also exhibits this property, we have run several simulations with the critical threshold parameter fixed to $L c , P , x − 1 = 1$ which, as discussed in Subsec. IV A, corresponds to a near-marginal state for $S 0 = 5$ [and f_d = 0]. In each run, the fueling rate has been set to a value in between 5 and 20. The resulting steady-state gradient profiles are shown in Fig. 5, where they can be seen to remain stiff over the central region between the sink and source, as expected, for $5 < S 0 < 10$⁠. At sufficiently larger fuelling rates, although the situation changes dramatically. For $S 0 = 20$⁠, the system is capable of driving a noticeable self-generated sheared poloidal flow (see Fig. 6) that inhibits turbulent transport near the source and sink, leading to a steepening of the averaged profile to compensate for the reduction of turbulence transport as clearly shown in Fig. 5. The average gradient actually steepens until it reaches values similar to those obtained in the absence of a critical gradient (i.e., for $L c , P , x − 1 = 0$⁠), as shown in Fig. 5, which implies that a sufficiently strong drive can develop self-generated flows capable of bringing the average profile back to a supermarginal state, where the gradient profile will be set by the source and the fluctuation levels, not the local threshold (⁠ $δ W P / ⟨ W P ⟩ t ≃ 20 %$ for $S 0 = 5$ and $δ W P / ⟨ W P ⟩ t ≃ 2.6 %$ for $S 0 = 20$⁠).

The last effect we have studied is whether the assumption (or not) of fast parallel equilibration changes the overall near-marginal transport dynamics of the model. We have carried out several simulations using $L c , P , x − 1 = 1 , L c , ⟨ P ⟩ y , x − 1 = 0 , L c , P , y − 1 = 0$⁠, and $S 0 = 5$⁠. These values correspond to a near-marginal steady-state for f_d = 0, as we showed in Secs. IV A and IV B. The fact that $L c , P , y − 1 = 0$ means that there is always free energy available to excite turbulence in the poloidal direction. The factor f_d has been varied between 0 and 1. The resulting steady-state gradient profiles are shown in Fig. 7. Curiously, the average gradients seem pretty similar at the two limiting values, f_d = 0 and f_d = 1, while the profile becomes closer to the supermarginal state (i.e., $L c , P , x − 1 = 0$⁠) when local and poloidally averaged threshold conditions have a similar importance (⁠ $f d = 0.4$⁠). The dynamics are however very different, as seen clearly when looking at the time (and poloidally) averaged poloidal flow (i.e., $⟨ v y ⟩ y , t$⁠), as shown in Fig. 8. Strong zonal flows play an important role in this change since they exist only in the sink and source regions for $f d ≪ 1$⁠, while the dominant threshold condition is the local one but is seen to extend more and more until they cover the whole domain for $f d ∼ 1$⁠, after the poloidally averaged threshold condition begins to dominate.

In this section, we will determine the fractional transport exponents that better captures the transport dynamics of each of the simulations discussed in Sec. IV A. In this way, we will be able to characterize quantitatively the nondiffusive nature of turbulent transport in the DTEM model as the dominance on transport dynamics is displaced from turbulent profile relaxation processes to zonal flows or vice versa. The methods described in Sec. III have been used for this task. Thus, the exponent H has been determined using the R/S analysis on the (x or y component of the) Lagrangian velocity time series of tracers advected by the turbulence.

The rescaled ranges (i.e., R/S) obtained for the tracer Lagrangian velocities of some simulations from the set in which the critical gradient threshold, $L c , P , x − 1$⁠, is varied are shown in Fig. 9 as a function of time lag. [As we mentioned earlier, for these simulations, $S 0 = 5$ and f_D = 0.] Two different scaling regions are apparent, separated at lag $τ Ω i ∼ 1$⁠. The region for $τ Ω i < 1$ exhibits a Hurst exponent close to 1, both for the x and y velocity components. This region is related to the autocorrelation of the time series with itself and tells us that single-eddy dynamics extend all the way up to timescales of the order of $τ Ω i ∼ 1$⁠. It is however irrelevant to determine the nature of long-term, long-distance transport dynamics. The second region, for $τ Ω i < 1$⁠, is the interesting one since, in a diffusive system, one should find $H ∼ 0.5$ [or, since we are using R/S, $H ∼ 0.55$] for timescales larger than single-eddy timescales. Instead, here we find that $H ∼ 0.8$ along the x direction for $L c , P , x − 1 = 0.8$ and $L c , P , x − 1 = 1$⁠, the cases with near-marginal gradient profiles, as we discussed in Sec. IV. This means that the nature of radial transport is superdiffusive, as one would expect from a near-marginal turbulent system in which transport is dominated by radial avalanches. On the other hand, transport is diffusive, or mildly subdiffusive, for $L c , P , x − 1 = 0$ and $L c , P , x − 1 = 0.4$⁠, since $H ∼ 0.49$⁠. This behavior is also consistent with the presence of supermarginal turbulence and the weakly radially sheared poloidal flows that we discussed in Sec. IV. Regarding the nature of poloidal transport for timescales $τ Ω i > 1$⁠, it is interesting to note that the Hurst exponent values are now reversed. $H ∼ 0.6$ (weakly superdiffusive) is obtained for $L c , P , x − 1 = 0$ and $L c , P , x − 1 = 0.4$⁠, while $H ∼ 0.5$ (diffusive, or weakly subdiffusive) for $L c , P , x − 1 = 0.8$ and $L c , P , x − 1 = 1$⁠. The persistence in the poloidal direction corresponds to the correlation induced by the self-generated flows.

In Self-Organized Criticality (SOC) systems, it is well-known that superdiffusion is only present for a certain range of scales in the R/S analysis, which is a consequence of the finite size of the system that can be traversed in a finite amount of time, above which correlations disappear.³⁰ In this periodic system, the system size is constrained between the source and the sink, and correlations tend to disappear once tracers travel more than several simulation boxes. The timescale in the range $1 < τ Ω i < 10$ corresponds to the average time for the tracers to travel from the source to the sink due to the relaxation events. Figure 9 also shows the appearance of this meso-range, which then corresponds to the near-marginal state as shown with in Fig. 3 and a reduction in the radially sheared, poloidal flows as shown in Fig. 4. In Sec. V B, when the fueling rate is varied, we will see that this superdiffusive range disappears when the tracers move across the source and the sink at a much faster rate, and the radial transport signature tends to $H ∼ 0.55$ for the largest S₀. Furthermore, the R/S analysis averages over multiple equal time segments, which means that there are multiple equal time segments of $1 < τ Ω i < 10$ that exhibit H > 0.55.

The flux-driven setup allows relaxation events, which are not present in fixed-gradient setups.^11,31 Although the “time-limited” subdiffusion and the transition to diffusion due to trapping and detrapping in coherent structures may persist, the dominant effect that determines the transition from superdiffusive to diffusive is associated with the lack of correlations due to finite-size effects. A flux-driven setup allows for the competition between relaxation events and self-generated sheared flows, which can be difficult to implement in a fixed-gradient setup. In this work, the R/S analysis averages over multiple equal time segments, which differs from the running diffusion coefficient that has been used to quantify the transition to diffusive transport in the other work.^11,31 In this work, the averaging of the time series per tracer is performed over the entire $t Ω i ∼ 500$⁠, which differs from the limited sampling time of $t Ω i ∼ 20$ associated with the previous observation.

We have collected the Hurst exponents along x and y obtained for all the simulations examined that varied $L c , P , x − 1$ in Fig. 10. The trends suggested by Fig. 9 become now much more apparent. With the uncertainty bars, the superdiffusive behavior along one direction appears to be correlated with the subdiffusive behavior along the other. This is expected due to the related transport characteristics observed in turbulence simulations with zonal-flows.⁷ In addition, it is clear that there is a minimum value of the critical threshold $L c , P , x − 1$ (close to 0.75) for near-marginal conditions to be established, in which avalanches dominate radial transport. Below this threshold value, zonal flow dynamics become dominant and reshape the nature of transport towards the type of subdiffusion characteristic of transport across zonal flows. The more subdiffusive, the stronger the radial shear of the poloidal flows is. A coarse parameter scan over $L c , P , x − 1$ for double the simulation resolution also shows a similar transition in the transport character (Fig. 10). However, the transition occurs at slightly a different parameter value due to both the increase in the source amplitude S₀ and a smaller dissipation scale. Since the self-generated poloidal flows are localized near the source and sink regions, with more spatial scales, the turbulent relaxation events are believed to disrupt the self-generated flows more effectively. This is then reflected by a shift to smaller values $L c , P , x − 1$⁠. As a result, a smaller value of the critical gradient is required to trigger a relaxation event.

We proceed next to determine the fractional exponent H for the set of simulations in which the source S₀ is varied. [As we mentioned earlier, for these simulations, $L c , P , x − 1 = 1$ and f_D = 0.] The rescaled ranges obtained for the tracer velocities of some selected simulations are shown in Fig. 11. Focusing again on the timescales of interest for transport (i.e., $τ Ω i > 1$⁠), we find that $H ∼ 0.8$ in the radial direction for $S 0 = 5 − 10$⁠, clear superdiffusive behavior. For $S 0 = 20$⁠, the radial Hurst exponent drops to $H ∼ 0.53$⁠, and thus, transport becomes diffusive or mildly subdiffusive. This abrupt change is consistent with what was described in Sec. IV B for the background gradient, which moved from submarginal to supermarginal conditions at $S 0 = 20$ due to the increasing excitation of larger zonal flows that could extend over the whole domain, away from the source and sinks. The behavior along the poloidal (i.e., y) direction also follows the trend previously described, with $H ∼ 0.51$ (diffusive or mildly subdiffusive) for $S 0 = 5 − 10$ and $H ∼ 0.58$ (mildly superdiffusive) for $S 0 = 20$⁠.

A more complete picture of these trends is shown in Fig. 12, where the values of H obtained for all simulations in this set are shown. Clearly, radial transport remains superdiffusive for $S 0 < 10 − 12$⁠. The fact that H varies very little in this range is a consequence (or a symptom) of the type of profile resiliency that is characteristic of near-marginal turbulent transport. As expected, transport in the perpendicular direction (y) is then diffusive or mildly subdiffusive. For $S 0 > 12$⁠, however, the situation begins to change, and for $S 0 > 15$⁠, it is completely reversed. Strong poloidal flows have become dominant, turbulence has become supermarginal, and radial transport becomes diffusive or mildly subdiffusive.

We describe next the values for the fractional exponent H obtained for the set of simulations in which f_D is varied from 0 to 1 [In these, $L c , ⟨ P ⟩ y , x − 1 = 0 , L c , P , x − 1 = 1$⁠; $L c , P , y − 1 = 0$⁠, and $S 0 = 5$] in order to explore the importance of parallel equilibration in the dynamics. Fig. 13 shows the rescaled ranges obtained for the tracer Lagrangian velocities in a few selected simulations. Regarding radial motion (i.e., along x), it is interesting to note that superdiffusive behavior is only found for f_d = 0 $( H ∼ 0.8 )$⁠, while strong subdiffusive behavior (H < 0.4) is found in all other cases shown. Consistently, with our previous findings, transport in the poloidal direction shows the opposite trend, with diffusive or mild subdiffusion for the case with f_d = 0 and superdiffusive motion for all others.

A more precise analysis can be done using Fig. 14, which collects the H values obtained for all runs made with varying f_D. The figure clearly shows that, as poloidal equilibration becomes sufficiently strong (i.e., for $f D > 0.2$ in this particular set), superdiffusive radial [poloidal] transport is abruptly replaced by subdiffusive [superdiffusive] transport. The more subdiffusive, the more dominant poloidal equilibration becomes, and the stronger and more radially extended, the poloidal flows become, as we discussed in Subsec. IV B.

It seems clear, from the results just described, that the interaction between self-consistently evolved background profiles and turbulence can yield a large variety of (radial) turbulent transport dynamics. On the one hand, in cases in which the background profiles manage to stay close to near-marginality at the steady-state and when flows remain small, transport behaves superdiffusively, as expected from the typical coherent, down- or up-the-gradient profile relaxation that is usually called “avalanche.” On the other hand, when the radial shear of the turbulence-driven poloidal flows manages to become large, they can not only decorrelate the coherent avalanches and turn the nature of radial transport more diffusive but also even make it behave subdiffusively. The physical mechanisms responsible for subdiffusion in this case have already been described somewhere and were due to the enhancement (or repression) of a certain sign of the parallel (to the magnetic field) fluctuating vorticity.

In the simulations presented in this paper, we have shown that these turbulent poloidal flows, being driven by the Reynolds stresses if sufficiently non-homogeneity and anisotropy exists in the turbulence,³² always appear first close to the source and sink region where gradients are large, particularly in the radial direction. In the presence of a constant source, however, these flows can be suppressed by increasing the critical threshold, $L c , P , x − 1$⁠, that reduces the level of fluctuations needed to balance the source, in spite of the larger inhomogeneity. This reduces the Reynolds stresses and the subsequent flows. At the same time, larger radial structures may develop (see Fig. 15), due to the shift to a larger growth rate. In this situation, it becomes easier that the steady-state of the system remains close to near-marginality, as we discussed in Sec. IV A, which makes coherent relaxations down-gradient (i.e., avalanches) more frequent. As a result, radial transport becomes superdiffusive, as confirmed by the variation of the H exponent found in our runs and shown in Fig. 10. The change in dynamics seems, however, to happen rather abruptly, which may point to a dynamical transition. In this context, we mean that a “dynamical transition” is an abrupt change in the transport signatures, which suggest a change in the dynamics rather than a divergent change in the sense of thermodynamical transition.

On the other hand, the opposite happens when the source S₀ is increased for a constant critical threshold, $L c , P , x − 1$⁠. The larger source increases the inhomogeneity in the near regions while, at the same time, results in larger levels of fluctuations in order to reach a balance between the larger source and the turbulent induced fluxes (Fig. 16). Consequently, it becomes easier for poloidal flows to be driven by the Reynolds stresses. If the initial steady-state was one of near-marginality, in which radial transport takes place mainly via radial avalanches, an increase in the external drive thus results in profiles becoming supermarginal and in (radial) transport becoming more diffusive or if the poloidal flows become sufficiently sheared in the radius, even subdiffusive. This is what the simulations shown in Sec. IV B and what the Hurst analysis confirmed as can be seen in Fig. 12. Again, the change in dynamics takes place rather suddenly, pointing to the possibility of a dynamical transition.

We can interpret the change in the radial dynamics further by looking at the evolution of the poloidal averaged radial flux $⟨ Γ P , x ⟩ y$ that can be calculated from $Γ P , x = P v x = − P ∂ y ϕ$⁠, where P > 0 always. With this convention, the sign of the flux refers to the flow direction. Fig. 16 shows the temporal evolution of these fluxes for increasing source S₀, from $S 0 = 5$ (frame (a)) to $S 0 = 20$ (frame (d)). In the near marginal state $S 0 = 5$⁠, diagonal features are apparent that reflect the ongoing coherent transport events (i.e., avalanches) that propagate down and up the gradient, which are reminiscent of the transport observed in critical-threshold sandpile models.¹⁹ The transition in the dynamics starts to be noticeable around $S 0 ∼ 10$ (Fig. 16(c)), where $∂ x P$ becomes sufficiently steep to generate sufficiently strong poloidal flows that can compete with the local relaxation process. In effect, the self-generated sheared flow introduces a decorrelating mechanism in the radial transport events induced by the local turbulence relaxation process. When the self-generated poloidal flows are enhanced by the larger source, correlated transport events become more scarce, which is reflected by a decrease in diagonal features as shown in Fig. 16(d). With the reduction in correlated flux events, $⟨ ∂ x P ⟩ y , t$ also becomes steeper on average as shown in Fig. 5.

Finally, we have also shown that the dynamics of zonal flow generation and saturation may be greatly impacted by the presence (or absence) of efficient parallel equilibration (along y, in our case). As a result, the overall radial transport dynamics is affected as well. A system whose steady-state profile is near-marginal in the absence of parallel equilibration (that is, when f_D = 0 in our model) can transit to supermarginal profiles due to the development of very large poloidal flows that are no longer mostly present at the source and sink locations, where the inhomogeneity is larger, but that may extend over the whole system. The reason here is that, by introducing parallel equilibration, the anisotropy has been increased enormously everywhere in the system. Indeed, the x and y directions, are now also distinguished by which is the dominant way in which they react to local gradients: by turbulent-induced fluxes in the x direction and by parallel equilibration along the y direction. As a result, larger radial structures develop (Fig. 2). The resulting larger anisotropy, in combination with the inhomogeneity introduced by the source, fuels the Reynolds stresses system-wide and generates the large poloidal flows whose action on the turbulence makes the profiles supermarginal, and radial transport becomes strongly subdiffusive, as was shown in Sec. IV C. The Hurst analysis also suggests that this change in nature is rather abrupt (happening at $f D ∼ 0.2$ in our case, as shown in Fig. 14), maybe reminiscent of a dynamical transition. Parallel equilibration also permits the formation of larger radial structures (see the third frame in Fig. 2), which set up a state in which down-gradient transfer from the source to the sink may be accomplished within a couple of eddies despite the self-generated poloidal flows across the entire profile, which also enhances subdiffusion.

In this paper, the interaction between self-consistently evolved profiles, turbulence, and turbulence-driven zonal flows has been explored in the context of a simplified drift-wave turbulence model. Although instances of non-diffusive behavior in plasma turbulence have been previously studied, the results obtained with our model are perhaps more compelling in the sense that they illustrate clearly how the resulting radial transport can be superdiffusive, diffusive, or subdiffusive depending on several factors. Namely, how strongly the system is driven (i.e., by changing S₀), the degree of effectiveness of turbulence to relax profiles (by changing $L c , P , x − 1$⁠, which sets the growth rate and how near profiles can be to marginal values), and the degree of competition offered by other transport mechanisms (parallel equilibration, by changing f_d = 0). The existence of this rich dynamical zoology is one of the main results of this paper, as well as the understanding gained on how one can transition from one regime to another in dynamical space by modifying the free parameters of the model.

The results presented here should also serve as a warning, once more, of the dangers of doing turbulent transport studies numerically by using fixed-gradient setups. The numerical evolution of turbulence in the presence of frozen background profiles often can produce diffusive transport. The interplay between turbulent relaxation and self-generated sheared poloidal flows, which form the basis for the transport explored in this model, is however absent unless a flux-driven setup is used. Indeed, most of the rich dynamics shown here were not present when running our simplified DTEM model without an equation for background profile evolution.

All of the above points to flux-driven, global gyrokinetic simulations as the more adequate tool to carry out similar studies in tokamak-relevant situations. The work presented here can serve as a guide of what should be expected, both from a physics and a methodological point of view. We feel that the study of the role of parallel equilibration in this context might be particularly relevant in tokamaks, particularly close to the plasma edge where collisional drag is a factor, or in the neighborhood of low-order rational surfaces, where equilibration along the field line may still yield significant on-surface variations.

This work was supported by U.S. DOE under Contract No. DE-FG02-04ER54741 with the University of Alaska Fairbanks and in part by a grant of HPC resources from the Arctic Region Supercomputing Center at the University of Alaska Fairbanks. This research was also sponsored in part by DGICT (Dirección General de Investigaciones Científicas y Tecnológicas) of Spain under Project No. ENE2015-68265.