2. Turbulent amplification process

In the combined SOC-turbulence scenario which we consider, a guiding role is attributed to the usual picture of eddies and eddy-induced transport associated with plasma instabilities. A paradigmatic framework with which to understand and explain the phenomenon from first principles is drift-wave turbulence, modelled by the HW equations (Hasegawa & Wakatani Reference Hasegawa and Wakatani1983), as direct numerical investigations have shown (Naulin, Nielsen & Rasmussen Reference Naulin, Nielsen and Rasmussen1999, Naulin et al. Reference Naulin, Garcia, Nielsen and Rasmussen2004; Basu et al. Reference Basu, Jessen, Naulin and Rasmussen2003). Nonlinearly, in a magnetic confinement geometry, the eddy-induced transport reduces the slope of the average profile where the vorticity is maximal, and, at the same time, steepens it in its nearby vicinity, thus increasing the instability in the next radial location. This displacement of the instability is an avalanche in that the step in the gradient moves radially outward, creating an unstable propagating frontFootnote ⁴ . However, because drift-wave turbulence is essentially two-dimensional, there exists an inverse cascade of the energy which is associated with the phenomena of eddy merging and the formation of large-scale coherent structures in the strongly turbulent flow. To this end, the propagation of the unstable front becomes a combined effect due to the next-radial generation of the off-spring eddies and their merging with the ever-growing mother eddy. One sees that turbulence will act to amplify the avalanches by fuelling them with more free energy via the inverse cascade. The process will stop when excess energy and particles are eventually let out through boundaries, thus reducing the slope of the average profile on the system-size scales.

We employ a simplification of the HW model (Hasegawa & Wakatani Reference Hasegawa and Wakatani1983) for two-dimensional electrostatic drift-wave turbulence driven by the resistive instability

Here the Poisson bracket $\{{\it\varphi},g\}=\hat{\boldsymbol{z}}\boldsymbol{\cdot }[\boldsymbol{{\rm\nabla}}{\it\varphi}\times \boldsymbol{{\rm\nabla}}g]$ is used to denote the nonlinear terms originating from advection with the $E\times B$ drift; $E=-\boldsymbol{{\rm\nabla}}{\it\varphi}$ is the intensity of the electrostatic field; ${\it\varphi}$ is the electrostatic potential; ${\rm\nabla}^{2}{\it\varphi}$ is the vorticity, that is, the curl of $E\times B$ ; and we have replaced the parallel derivatives with an effective parallel wavelength. A damping term ${\it\kappa}{\rm\nabla}^{4}{\it\varphi}$ is added to (2.2) to dissipate vorticity on the small scales. It is assumed that only one mode parallel to the homogeneous outer magnetic field $\boldsymbol{B}=B_{0}\hat{\boldsymbol{z}}$ is excited, so that all differential operators work in the $(x,y)$ -plane. We may identify the $x$ -coordinate with the radial direction and the $y$ -coordinate with the poloidal direction, as a portion of the plasma-edge in a toroidal geometry is represented. Further ${\it\delta}=1/k_{\Vert }^{2}L_{\Vert }^{2}$ is the adiabaticity parameter and characterizes the deviation between the potential and the density fluctuations in the HW model. In a basic theory of drift-waves, it is shown that this deviation leads to an instability with a maximum linear growth rate ${\it\gamma}_{L}\approx {\it\delta}/8$ at $k={\rm\nabla}^{-1}\approx 1$ . It is therefore convenient to think of ${\it\delta}$ as a driving rate for the turbulence. The inverse adiabaticity parameter, $1/{\it\delta}$ , characterizes the coupling strength (in the complexity sense of the wording) between the two interacting degrees-of-freedom represented by, respectively, the ${\tilde{n}}$ and ${\it\varphi}$ dependencies. It should be noted that ${\it\delta}$ absorbs via the scale length $L_{\Vert }=(L_{n}T_{e}/m_{e}c_{s}{\it\nu}_{ei})^{1/2}$ the parameters of parallel dynamics; it also contains through the electron–ion collisional frequency ${\it\nu}_{ei}$ the parallel resistivity. Here, $L_{n}$ is the scale length of the perpendicular background density gradient: $L_{n}^{-1}=|\boldsymbol{{\rm\nabla}}n/n_{0}|$ , where $n=n_{0}(x)+{\tilde{n}}(x,y,t)$ . The normalization scales are $c_{s}/{\it\Omega}_{i}$ for lengths perpendicular to the magnetic field and $L_{n}/c_{s}$ for the time, where ${\it\Omega}_{i}$ is the ion cyclotron frequency and $c_{s}=\sqrt{T_{e}/m_{i}}$ is the sound speed. Although the mathematical structure of the HW equations (2.1) and (2.2) is clear and simple, their solutions for turbulent flows are immensely complex and can only be obtained approximately through numerical computation. This absence of simple solutions is due to the nonlinear nature of the equations, giving rise to mode coupling, build-up of correlations and the formation of coherent vortical structures in the strongly turbulent flow.

We are now in a position to formulate a simple criterion for the onset of avalanching dynamics and its amplification by inverse cascade of the energy. This criterion shall use the idea of separation between wave-like motions and nonlinear structures. In drift-wave turbulence, the so-called Rhines length, ${\it\lambda}_{Rh}$ , determines the upper bound on the size of vortical structures in the flow. The Rhines length, originally introduced in geophysical fluid turbulence (Rhines Reference Rhines1975), and later applied to drift-wave turbulence (Naulin Reference Naulin2002), designates the spatial scale separating vortex motion from drift-wavelike motion. For the HW system, it is obtained as ${\it\lambda}_{Rh}\propto \sqrt{u_{E\times B}}$ , leading to a characteristic turnover time ${\it\tau}_{turn}\sim {\it\lambda}_{Rh}/u_{E\times B}\propto 1/\sqrt{u_{E\times B}}$ (termed the Rhines time). The avalanching behaviour occurs when the eddy turnover time, ${\it\tau}_{turn}$ , is small compared with the instability growth time. For the purpose of formal ordering, we require ${\it\tau}_{turn}\ll {\it\gamma}_{L}^{-1}$ , where ${\it\gamma}_{L}$ is the linear growth rate. Using here that the linear growth rate is proportional to the non-adiabaticity of the fluctuations, i.e. ${\it\gamma}_{L}\propto {\it\delta}$ , we can order, up to numerical coefficients,

One sees that the turbulence driving rate, represented by the parameter ${\it\delta}$ , must be small compared with the inverse Rhines time in the system. The time scale separation in (2.3) is favoured in the regime of strong coupling or near adiabaticity, on the one hand, and in the presence of large electrostatic fluctuations or fast $E\times B$ drifts, on the other hand. In the latter case, the restrictive assumptions of smallness of the fluctuations underlying the HW model are invalidated. Furthermore, since the instability drive of the turbulence is controlled by electron–ion collisions in the parallel electron motion via ${\it\delta}\propto {\it\nu}_{ei}/k_{\Vert }^{2}$ , the condition for time scale separation in (2.3) implies that the particle cross-field transportation time is short compared with the characteristic resistive time. Thus, the avalanches tend to dissipate their content mainly upon reaching the boundaries, and merely redistribute the particles and free energy across the system otherwise. These conditions of slow driving and time scale separation, combined with the stabilizing role of boundary dissipation, constitute a typical set-up for dynamical systems exhibiting SOC (Aschwanden Reference Aschwanden and Aschwanden2013; Aschwanden et al. Reference Aschwanden, Crosby, Dimitropoulou, Georgoulis, Hergarten, Jensen, McAteer, Milovanov, Mineshige and Morales2014). At this point, the propagation of unstable fronts due to the processes of eddy merging and interactions acquires the typical signatures of the avalanching dynamics of the SOC type. Nonlinearly, the absorbing boundary at the plasma edge provides a feedback of the particle and energy loss processes on the dynamical state of bulk plasma; where the system of gradients pumping the drift-waves self-adjusts to generate self-organization in a marginally stable state. Thus, SOC attracts nonlinear feedback dynamics (Sornette Reference Sornette1992). Moreover, the condition for time scale separation in (2.3) suggests that conservation of the transported quantity is indeed necessary to obtain criticality (e.g. Zhang Reference Zhang1989). It is understood that dissipative boundaries are perceived as an essential key element of the model as they introduce a feedback dynamics generating SOC.

3. The Galton–Watson chain process near extinction

Here we refine the avalanching model above and employ the idea of stimulated vortex formation (Milovanov & Rasmussen Reference Milovanov and Rasmussen2014). This means that the original process of vortex formation occurs through the stimulated birth of one vortex due to a nearby recently formed vortex, and that the process, which is chain reaction-like, is self-controlling and reduces itself via nonlinear feedback dynamics to its marginal instability threshold. We think of the turbulence as starting from a uniform fluctuation background within which there is some small probability of spontaneous vortex formation. Then the birth of one vortex will perturb the unstable background nearby such that one or more vortices will appear next to the mother vortex. In two-dimensions, the energy cascading from smaller to larger scales will cause the mother and the off-spring vortices to merge, and the process will repeat itself. This process of stimulated vortex formationFootnote ⁵ will generate an ever-amplifying instability front until the transferred energy is let out through the boundaries. The amplification of the vortical motions finds its energy reservoir in the non-adiabatic properties of the fluctuations, characterized by the parameter ${\it\delta}$ .

We assume that the original process of vortex formation occurs near marginality and is self-reproducing. Then, it can be modelled as a variant of the Galton–Watson chain process near extinctionFootnote ⁶ . This process will be self-similar in $2+1$ dimensions, with time interpreted as the preferred dimension (see figure 1). Note that the ‘arrow of time’ looks in the direction of the energy transfer (the direction of amplification). Hence the time coordinate plays a very special role in the model in that it introduces in a parametric form a subordination to the hierarchy of vortices. The spatio-temporal character of the Galton–Watson process matches with the implication of SOC, where the space correlations in hierarchic geometry act to attract the nonlinear feedback dynamics (Sornette Reference Sornette1992; Milovanov Reference Milovanov and Aschwanden2013).

Figure 1.

Stimulated vortex formation as a Galton–Watson chain process in $2+1$ embedding dimensions, with time, $t$ , interpreted as the preferred dimension. An artist’s view, using a continuum background to accentuate the spatio-temporal character of the dynamics behind the phenomena of SOC. The ‘arrow of time’ looks in the direction of energy transfer. Adapted from Milovanov & Rasmussen (Reference Milovanov and Rasmussen2014).

A hierarchic process being the Galton–Watson chain process near extinction implies that its size probability mass function is a pure power-law: $\hat{{\it\chi}}({\it\lambda})\sim {\it\lambda}^{-{\it\tau}+1}$ . The exponent ${\it\tau}$ , which is also introduced in major SOC theories (e.g. Tang & Bak Reference Tang and Bak1988; Zhang Reference Zhang1989), is obtained here using a binomial distribution and the Catalan numbers (Corral & Font-Clos Reference Corral, Font-Clos and Aschwanden2013). The end result is ${\it\tau}=5/2$ and is mathematically exact. This value is robust in that it does not depend on the specific version of the Galton–Watson process employed. In the above, we used hat to emphasize that the probability distribution $\hat{{\it\chi}}({\it\lambda})$ is defined in $(2+1)$ -dimensional space–time. In the limit of long times (small wave vectors), however, one can factorize the dependencies in $\hat{{\it\chi}}({\it\lambda})$ to obtain ${\it\chi}({\it\lambda})\sim \hat{{\it\chi}}({\it\lambda})/{\it\lambda}\propto {\it\lambda}^{-{\it\tau}}$ , where ${\it\chi}({\it\lambda})$ is the usual eddy size distribution in the real space. Assuming that there is a characteristic flow velocity, the $E\times B$ velocity, the distribution of eddy sizes can be translated into a heavy-tailed flux distribution, ${\it\chi}({\it\Gamma})\sim {\it\Gamma}^{-{\it\tau}}$ , where ${\it\Gamma}=\int \,\text{d}{\it\sigma}[{\tilde{n}}u_{E\times B}]/\!\int \,\text{d}{\it\sigma}$ is the average particle flux distribution and is obtained as a flux-surface integral of the convected density. This, together with the fact that the turbulent perpendicular fluxes are determined by the trapping and subsequent displacements of particles with vortical motions, leads to the inverse-power Pareto distribution of jump lengths

where the exponent of the power-law is $-5/2$ exactly.

4. Lévy-fractional diffusion equation

The inverse power-law dependence in (3.1) is characteristic of Lévy stable motion processes often referred to as Lévy flights. The Lévy flights represent perhaps the simplest toy model for non-local transport, and in the absence of external force fields, correspond to an asymptotic ( $t\rightarrow +\infty$ ) diffusion equation (e.g. Metzler & Klafter Reference Metzler and Klafter2000 for review)

where $1<{\it\mu}<2$ is the Lévy index and

is the Riesz operator in infinite space. The dependence in front of the second-order space derivative has the sense of a numerical normalization parameter. Setting an absorbing boundary at $x=0$ , with an initial condition defined along the positive semi-axis $x>0$ , one requires that the solution necessarily has to vanish on the negative semi-axis. Then the non-local integration from $-\infty$ to $+\infty$ in (4.2) is reduced to the positive semi-axis, motivating

from which the first passage time density for Lévy flights can be deduced (Chechkin et al. Reference Chechkin, Metzler, Gonchar, Klafter and Tanatarov2003; Metzler & Klafter Reference Metzler and Klafter2004). Moreover, the model in (4.3) observes the Sparre Andersen universality (Sparre Andersen Reference Sparre Andersen1953) that the first passage time density decays as ${\sim}t^{-3/2}$ after $t$ time steps ( $t\rightarrow +\infty$ ). This behaviour has been confirmed through numerical simulation in Chechkin et al. (Reference Chechkin, Metzler, Gonchar, Klafter and Tanatarov2003). More so, starting from Lévy-fractional diffusion equations in (4.1) and (4.3), one infers for large ${\rm\Delta}x$ an asymptotic jump length p.d.f. ${\it\chi}({\rm\Delta}x)\sim |{\rm\Delta}x|^{-{\it\mu}-1}$ . This is equivalent with the distribution in (3.1), if

Using ${\it\tau}=5/2$ , one obtains ${\it\mu}=3/2$ . We associate this value with the phenomena of non-local transport generated intrinsically via the mechanisms of SOC-turbulence coupling through the inverse cascade. Non-local implies that the transport occurs on multiple scales and that the dispersion is so broad that the second moments diverge at long jumpsFootnote ⁷ . This situation of diverging second (and higher) moments is different from classic reductionism of Fickian diffusion, which is based on the assumption that the characteristic length exists and is finite. Indeed, from (3.1) one sees that $\int |{\rm\Delta}x|^{2}{\it\chi}({\rm\Delta}x)\,\text{d}{\rm\Delta}x\propto |{\rm\Delta}x|^{1/2}\rightarrow +\infty$ for ${\rm\Delta}x\rightarrow \infty$ ; whereas for a Gaussian ${\it\chi}({\rm\Delta}x)$ it converges at a fast pace.

That drift-wave turbulence in magnetic confinement devices can generate motion processes of the Lévy type has been discussed based on numerical simulations by several authors (del-Castillo-Negrete et al. Reference del-Castillo-Negrete, Carreras and Lynch2004; del-Castillo-Negrete Reference del-Castillo-Negrete2006), however, the parameter range of non-local transport has not been established. Nor has the fundamental significance of ${\it\mu}$ been understood, which remained in those simulations simply a numerical fitting index. Here, we find that non-local behaviour occurs in the regime of strong nonlinearity and time scale separation, (2.3), actually stretching values of parameters to the margins of validity of the HW system. Also strong nonlinearity implies fast $E\times B$ drifts, hence appreciable departures from the state of marginal stability. So the paradigmatic set-up based on the smallness of fluctuations and perturbation-theory analysis proves to be relaxed as soon as Lévy processes are concerned, paving the way to non-perturbative theories of anomalous transport (Milovanov & Iomin Reference Milovanov and Iomin2014, Reference Milovanov and Iomin2015). Another important point worth noting is that the HW system alone is not enough to justify non-local transport models, if only due to the restrictive assumptions of the smallness of fluctuations as well as the peculiar physical implications behind the processes of stimulated vortex formation and turbulent amplification of the avalanches. Indeed, the power-law p.d.f. in (3.1) uses as an essential ingredient the fact that there is an amplification process taking place, which we associate with the back-reaction of inverse cascade of the energy on SOC. It is in fact the combined effect of SOC-like and turbulence-like phenomena that produce, via complexity coupling, a non-local diffusion in the sense of Lévy-fractional equations (4.1) and (4.3). In this interpretation, the fractional exponent ${\it\mu}=3/2$ collects in one distinguished value the strong nonlinearity of the coupled SOC-turbulence system. It leads to a non-trivial description in which the kinetic equations (4.1) and (4.3) are formally linear, with all nonlinearities already absorbed by the fractional Laplacian in (4.2).

5. Spectral energy density

Based on the above analysis, we might suggest that the mechanism of turbulent amplification of the SOC avalanches is responsible for the occurrence of large-amplitude transport events in plasma confinement, associable with large intermittent bursts of transport. One question of practical importance is concerned with the precursors of this behaviour. Because of amplification, we expect a steeper drop-off in the energy spectrum of the coupled SOC-turbulence system as compared to the inertial range of the fluid (drift-wave) turbulence. Using here the fact that the divergence of the $E\times B$ drift is zero, which keeps the flow incompressible, one finds after many merging events that the variance

will be proportional to the square of the vortex size, leading to $\text{Var}\,u_{E\times B}\,({\it\lambda})\propto {\it\lambda}^{2}$ . When translated into wave vectors, ${\it\lambda}\sim 1/k$ , this becomes $\text{Var}\,u_{E\times B}\,(k)\propto k^{-2}$ . Employing for the spectral energy density, $S(k)$ ,

one readily obtains $S(k)\sim k^{-3}$ . Indeed, the spectrum of measured fluctuations in tokamak plasmas involves a $k^{-3}$ power-law subrange already discussed in Reuss & Misguich (Reference Reuss and Misguich1996). By contrast, assuming that the energy cascade rate does not depend on $k$ , one arrives at the scaling $\text{Var}\,u_{E\times B}\,(k)\propto k^{-2/3}$ , from which the familiar Kolmogorov spectrum $S(k)\sim k^{-5/3}$ is inferred. In this connection, we should stress that the avalanching transport is triggered by the explicit radial dependence in the profiles, and, when account is taken for the inverse cascade, by boundary feedbacks, so the assumptions of constant energy transfer and of infiniteness of the system, resulting in the fluid-like $-5/3$ behaviour, do not really apply here.

6. Summary and concluding remarks

We first list what we believe are the key features of the turbulent amplification process illustrated above. Clearly, transport is multi-scale and, as a general situation, non-local (in the sense of Lévy–Gnedenko generalized central limit theorem). Non-locality occurs naturally in self-consistent strong turbulence through the very specific process of stimulated vortex formation in the presence of the inverse cascade of the energy in reciprocal space. But it only comes on the scene under the condition that nonlinearity is strong and is indeed above a certain critical level defined as a competition between the Rhines time ${\it\tau}_{Rh}\propto 1/\sqrt{u_{E\times B}}$ and the instability growth time. Coupling between the unstable modes triggers the formation of coherent structures, which, under the condition for time scale separation, occur in Galton–Watson chains via the processes of the stimulated birth. Chain processes near extinction attract the dynamics of the vortexes in much the same way as standard SOC attracts the dynamics of the wave-like fluctuations to the state of marginal stability. The processes of amplification taking place imply enhanced probabilities of big events and may also favour events of anomalously large size beyond the range of validity of pure SOC schemes.

In the above discussion it was argued that a pure SOC state does not occur when the dominant transport mechanism in the bulk plasma involves other than the drift-wavelike patterns, the formation of coherent structures. So the dynamical system with coupling between the modes will evolve into a non-equilibrium, hybrid state in which the SOC and the coherent components operate in concert to mediate the large levels of transport to the boundaries. Similar hybrid SOC-coherent states have also been proposed for the plasma sheet of the Earth’s magnetotail (Milovanov Reference Milovanov and Aschwanden2013, pp. 161–162). The remarkable universality of such phenomena in plasma turbulence suggests they possibly have the same nonlinear origin. When a saturated turbulence state is approached, the initial fine-scale signatures of the fluctuations (both potential and the plasma density) will nonlinearly be washed out by amplification. When this occurs, the dynamics are dominated by the bulk-average nonlinearities, and not anymore by multi-scale features of the fluctuations in the HW regime, permitting appreciable departures from the state of marginal stability. Analysis of these strongly nonlinear regimes with feedback requires global models accounting for a self-consistent evolution of the background density profiles alongside the phenomena of mode coupling and the spectral energy transfer.

We should stress that the onset of non-local transport is limited to strong nonlinearities permitting both the formation of coherent structures (via the phenomena of mode coupling and the build-up of correlations) and the Galton–Watson chain process of stimulated vortex birth. It is these conditions, represented in a synthesized form by (2.3), which justify the introduction of the Riesz derivative into the transport model in (4.3) above. The physics behind this equation involves the idea of turbulent amplification of SOC by the inverse cascade of the energy, and in this sense, represents SOC-like and fluid turbulence-like features on essentially an equal footing. By contrast, pure SOC schemes lead to a different description as they correspond to kinetic equations with a fractional time, rather than space, derivative, owing to the next-neighbour character of the interactions (Milovanov Reference Milovanov2009, Reference Milovanov2011; Milovanov & Iomin Reference Milovanov and Iomin2012).

The exponent of fractional differintegration ${\it\mu}$ in (4.1) determines the characteristic anomalous scaling between the lifetime of the avalanches and their size, i.e. $t\sim |x|^{{\it\mu}}$ . Often such scaling relations are found phenomenologically in fluctuations of the low confinement mode plasma, where they are interpreted in terms of SOC (e.g. Sharma et al. Reference Sharma, Aschwanden, Crosby, Klimas, Milovanov, Morales, Sanchez and Uritsky2015). Indeed, the scaling $t\sim |x|^{{\it\mu}}$ is achieved in major SOC theories using the assumption of locality and the details of the lattice redistribution rule. In those theories, the exponent ${\it\mu}$ , which is model dependent, determines the universality class of SOC (Maslov & Zhang Reference Maslov and Zhang1996; in their notation, ${\it\mu}=z$ ). There is no indication that our mixed SOC-turbulence model falls within the known universality classes. Employing the concept of the Galton–Watson chain process near extinction, and the premise of stimulated vortex formation, we find ${\it\mu}=3/2$ exactly. The model observes the Sparre Andersen universality in the presence of an absorbing boundary at the plasma edge. We consider the Sparre Andersen universality as a characteristic property of the avalanching system with amplification.

The idea that self-consistent strong turbulence may amplify SOC via coupling to the inverse cascade is theoretically very general and is not restricted to (although greatly motivated by) the realm of fusion research. Here we promote this idea as a proxy to understand the catastrophic events in complex systems: severe magnetospheric storms, climate disruptions, huge solar flares and outstanding coronal mass ejections (Aschwanden et al. Reference Aschwanden, Crosby, Dimitropoulou, Georgoulis, Hergarten, Jensen, McAteer, Milovanov, Mineshige and Morales2014). It has been proposed recently that the back-reaction of the inverse cascade of the energy on SOC, mediated by drift-waves, poses a significant threat to big fusion devices such as future power plants (Milovanov & Rasmussen Reference Milovanov and Rasmussen2014), and that some catastrophic events may have clear precursors in the statistical datasets. The topic is of interest from a fundamental scientific perspective and for the practical realization of diagnostic tools to monitor, forecast and prevent large-amplitude bursts of transport in magnetic fusion reactors.

It is understood that forecasting big events is not at all trivial, partly due to their respective rareness and hence the lack of statistically significant data, and partly due to the problem of detection of the ‘dragon-kings’ of the statistics (Sornette & Ouillon Reference Sornette and Ouillon2012). The dragon-kings are catastrophic events of anomalously large size, lying off the extrapolation metrics in the corresponding ‘simpler’ systems (e.g. small-size tokamaks), and their correct identification and detection through the datasets remains a very delicate and difficult field, if only due to the scarcity of data as well as the extraordinary important implications with respect to hazard assessment, risk control and predictability. In many ways, the phenomenon of the dragon-king echoes Anderson’s ‘More and different’ (Reference Anderson2011) and in the context of controlled fusion, brings up, rephrasing de Gennes:

The idea (‘to put the sun into a box’) is pretty. The challenge is that the sun ought not to remain in a box.

That non-local features can occur naturally via self-organization could be demonstrated based on self-tuning phase transitions, using the idea of competing order parameters and the formalism of the fractional Ginzburg–Landau equation (Milovanov & Rasmussen Reference Milovanov and Rasmussen2005; Milovanov Reference Milovanov and Aschwanden2013). This may indeed prove to be a very promising (and quite novel using the notion of Galton–Watson chain process near extinction) angle of attack to characterize the origin of non-local transport in magnetically confined plasma systems. When seen from the perspective of time scale separation, the condition for strong nonlinearity favours the idea that strong (vortical dominated) drift-wave turbulence can acquire features that are otherwise associable with SOC. Indeed, the patterning adheres to a situation in which SOC-like and turbulence-like phenomena act in concert to regulate particles and heat transport to the boundaries, and are not separable statistically. Other than fusion, these phenomena of SOC-turbulence coupling might be proposed for geophysical flows (using as appropriate the notion of the Rhines length), where they shall be responsible for outstanding perturbations beyond the expected weather and climate patterns. Research on these topics using statistical data sets from climate dynamics might be strongly advocated.