2.1. Fixed gradient Hasegawa–Wakatani model

We consider the Hasegawa–Wakatani model in a 2-D slab $(\boldsymbol{e}_{x},\boldsymbol{e}_{y})$ of size $L_{x}\times L_{y}$ orthogonal to a uniform and constant magnetic field $\boldsymbol{B}=B\boldsymbol{e}_{z}$ , where $x$ and $y$ are respectively the radial and poloidal directions normalised to the sound Larmor radius $\rho _{s}$ . We first assume that the fluctuations are periodic in both axes.

The modified Hasegawa–Wakatani equations are (Hasegawa & Wakatani Reference Hasegawa and Wakatani1983; Numata et al. Reference Numata, Ball and Dewar2007; Guillon & Gürcan, Reference Guillon and Gürcan2025)

(2.1a) \begin{align} \partial _{t}\varOmega +\left [\phi ,\varOmega \right ] & =C(\widetilde {\phi }-\widetilde {n}),\\[-10pt]\nonumber \end{align}

(2.1b) \begin{align} \partial _{t}n+\left [\phi ,n\right ]+\kappa \partial _{y}\phi & =C(\widetilde {\phi }-\widetilde {n}), \end{align}

where $\phi$ is the electrostatic potential, $n$ the density perturbation and $\varOmega ={\nabla} _{\perp }^{2}\phi$ corresponds to the vorticity. The Poisson bracket $[\phi ,A]=\partial _{x}\phi \partial _{y}A-\partial _{y}\phi \partial _{x}A$ represents the advection of $A$ by the $E\times B$ velocity $\boldsymbol{v}_{E\times B}=\boldsymbol{e}_{z}\times \boldsymbol{\nabla }\phi$ . The adiabaticity parameter $C= {k_{||}^{2}\rho _{s}^2B}/{en_{0}\eta _{e,||}}$ , with $e$ the unit charge, $k_{||}$ the parallel wavenumber of fluctuations (along the magnetic field lines), $\eta _{e,||}$ the parallel resistivity and the background density gradient $\kappa \equiv - ({1}/{n_{r}}) ({\text{d}n_{r}}/{\text{d}x})$ , with $n_{r}$ being the background density profile normalised to a constant density reference $n_{0}$ , are constant parameters of the model. Even though here, we consider the inviscid case for simplicity, dissipation terms are generally introduced in both equations for numerical studies. Note that we have also decomposed the perturbations $A$ into their zonal $\overline {A}$ (averaged over the poloidal direction) and non-zonal $\widetilde {A}$ components, with a poloidal average defined as

(2.2) \begin{equation} \overline {A}(x,t)=\langle A\rangle _{y}\equiv \frac {1}{L_{y}}\int _{0}^{L_{y}}A\,\text{d}y,\ \langle \widetilde {A}\rangle _{y}=0\,. \end{equation}

In fact, (2.1b ) can also be written for the total density $N(x,y,t)=n_{r}(x,t)+\widetilde {n}(x,y,t)$ , normalised here to the constant reference density $n_{0}$ , so that

(2.3) \begin{equation} \partial _{t}N+\left [\phi ,N\right ]=C(\tilde {\phi }-\tilde {n}), \end{equation}

where

(2.4) \begin{equation} n_{r}=n_{lin}+\overline {n} \end{equation}

is the radial density profile, which combines the zonal density fluctuations with the linear (i.e. fixed gradient) background density profile as

(2.5) \begin{equation} n_{lin}(x,t)=-\kappa (x-L_{x})+n_{r}(x=L_{x})\,. \end{equation}

The Poisson bracket for $n_{lin}$ indeed yields $[\phi ,n_{lin}]=-\partial _{x}n_{lin}\partial _{y}\phi =\kappa \partial _{y}\phi$ .

In a fixed gradient formulation, since $\kappa$ is a constant, the background density profile $n_{lin}$ is prescribed and does not change in time. One may argue that this corresponds to imposing a complicated form of source terms to the density profile, which compensates the turbulent particle flux generated by the underlying instability (Garbet & Waltz Reference Garbet and Waltz1998), and that the remaining zonal density is basically a radially periodic perturbation, which cannot change the mean gradient.

2.2. Flux-driven model

To develop the flux-driven formulation, let us first consider $n_{r}$ to be an arbitrary density profile subject to a transport equation due to the turbulent particle flux.

2.2.1. Decomposition of the radial density profile

We can still decompose the density profile into a ‘linear’ part $n_{lin}$ and a radially periodic perturbation part $\overline {n}$ , which is just the zonal part of the density profile

(2.6) \begin{equation} n_{r}(x,t)=n_{lin}(x,t)+\overline {n}(x,t), \end{equation}

where

(2.7) \begin{equation} n_{lin}(x,t)\equiv -\kappa (t)(x-L_{x})+n_{r}(x=L_{x},t)\, \end{equation}

is no longer a constant in time, so that

(2.8) \begin{equation} \kappa (t)\equiv -\frac {n_{r}(x=L_{x},t)-n_{r}(x=0,t)}{L_{x}}\, \end{equation}

is a function of time, defined for the radial profile at a given time, and corresponds to the mean density gradient in the radial direction. The decomposition is illustrated in figure 1.

Figure 1. Decomposition of a given radial profile $n_{r}$ (blue) into its linear component $n_{lin}$ (dashed green) and its zonal fluctuating component $\overline {n}$ (dash-dotted red). Notice that using such a decomposition, we find that the zonal profile is periodic at the edges, while its derivative $\partial _{x}\overline {n}$ is not.

This decomposition is important because it enables us to obtain a radially periodic perturbated part, i.e. $\overline {n}$ , from which we will compute the fast Fourier transform, while the linear part will correspond to the constant in space $\kappa (t)$ in (2.1b ), which is now time dependent. Furthermore, this decomposition may help us project what we know about the linear, zonal and transport properties of the classic gradient-driven system into the flux-driven version, where the linear parameter $\kappa$ will now evolve in time in response to the turbulent particle transport. We know for example from fixed gradient studies of the Hasegawa–Wakatani system (Numata et al. Reference Numata, Ball and Dewar2007; Grander et al. Reference Grander, Locker and Kendl2024; Guillon & Gürcan Reference Guillon and Gürcan2025) that for $C/\kappa \lesssim 0.1$ , we get 2-D-like turbulence dominated by eddies, while for $C/\kappa \gtrsim 0.1$ , zonal flows dominate. We can argue that this observation should apply equally using a time evolving local mean gradient in a flux-driven system.

Notice that we can also define the zonal density to have a zero mean value:

(2.9) \begin{equation} \overline {n}\equiv n_{r}-n_{lin}-\left \langle n_{r}-n_{lin}\right \rangle _{x}\!, \end{equation}

where $\langle \boldsymbol{\cdot }\rangle _{x}$ denotes averaging along the $x$ axis. Thus, we have

(2.10) \begin{equation} n_{r}=n_{lin}+\overline {n}+n_{m}, \end{equation}

with $n_{m}\equiv \langle n_{r}-n_{lin}\rangle _{x}$ , so that $\overline {n}$ is really a zero mean, radially periodic perturbation associated with the density profile, as it would be the case in a fixed gradient model.

2.2.2. The equations

We can then write the time derivative of the total density $N(x,y,t)=n_{r}(x,t)+\widetilde {n}(x,y,t)$ using (2.3) as

(2.11) \begin{align} \partial _{t}N+\left [\phi ,N\right ] & =\partial _{t}n_{r}+\partial _{t}\widetilde {n}+\left (\kappa (t)-\partial _{x}\overline {n}\right )\partial _{y}\phi +\left [\phi ,\widetilde {n}\right ]=C(\widetilde {\phi }-\widetilde {n}), \end{align}

where we decomposed the Poisson bracket of the total density. Averaging this along $y$ , we get the evolution equation for the radial density profile:

(2.12) \begin{equation} \partial _{t}n_{r}+\partial _{x}\varGamma _{n}=0, \end{equation}

with $\varGamma _{n}$ the radial particle flux that is averaged along $y$ , defined as

(2.13) \begin{equation} \varGamma _{n}(x,t)\equiv \langle \widetilde {n}\widetilde {v}_{x}\rangle _{y}=-\langle \widetilde {n}\partial _{y}\widetilde {\phi }\rangle _{y}\,. \end{equation}

Note that density source terms can be introduced in (2.12). Subtracting the averaged (2.12) from the complete (2.11), we get the time evolution of the non-zonal density fluctuations

(2.14) \begin{equation} \partial _{t}\widetilde {n}+\left (\kappa (t)-\partial _{x}\overline {n}\right )\partial _{y}\phi +[\phi ,\widetilde {n}]-\langle [\widetilde {\phi },\widetilde {n}]\rangle _{y}=C(\widetilde {\phi }-\widetilde {n}), \end{equation}

so that we can write the complete flux-driven model as

(2.15a) \begin{align} \partial _{t}\varOmega +\left [\phi ,\varOmega \right ] & =C(\widetilde {\phi }-\widetilde {n})+\nu {\nabla} ^{2}\widetilde {\varOmega }, \end{align}

(2.15b) \begin{align} \partial _{t}\widetilde {n}+\left (\kappa (t)-\partial _{x}\overline {n}\right )\partial _{y}\widetilde {\phi }+\partial _{x}\overline {\phi }\partial _{y}\widetilde {n}+[\widetilde {\phi },\widetilde {n}]-\langle [\widetilde {\phi },\widetilde {n}]\rangle _{y} & =C(\widetilde {\phi }-\widetilde {n})+D{\nabla} ^{2}\widetilde {n}, \end{align}

(2.15c) \begin{align} \partial _{t}n_{r}+\partial _{x}\varGamma _{n} & =S_{n}(x,t)+D_{0}\partial _{x}^{2}n_{r}, \end{align}

where $\kappa$ and $\varGamma _{n}$ are given by (2.8) and (2.13), respectively. Note that viscosity and particle diffusion coefficients $\nu$ and $D$ , which act on vorticity and density fluctuations, respectively, have been introduced, but they could as well be replaced by hyperdiffusivity terms. An explicit diffusion term $D_{0}\partial _{x}^{2}n_{r}$ , which represents neoclassical diffusion on the radial profile, has also been added, along with a source term $S_{n}(x,t)$ , which represents particle injection.

We could apply the same decomposition procedure to a given radial profile $u_{r}$ of the $E\times B$ velocity, which would allow us to introduce the effect of large-scale mean sheared flows, generated by the radial electric field $E_{r}$ , which plays a fundamental role in turbulence suppression and pedestal formation in the L–H transition (Biglari, Diamond & Terry Reference Biglari, Diamond and Terry1990; Hinton & Staebler Reference Hinton and Staebler1993a ). However, this is technically more complicated because the vorticity (2.15a ) features the first and the second derivatives of the radial profile of the electrostatic potential $\phi _{r}$ (respectively the velocity $u_{r}=\partial _{x}\phi _{r}$ and the vorticity $\varOmega _{r}=\partial _{x}^{2}\phi _{r}$ radial profiles), and it is not clear which decomposition between a linear and zonal profile would be the better adapted.Footnote ¹ We leave this discussion to a future publication where an external $E_{r}$ shear acts as an important control parameter.

Having presented the theoretical model, we now turn our attention to the numerical implementation of the FDHW system, using a pseudo-spectral method with penalisation.

3.1. Buffer zone

In the flux-driven formulation, in addition to all the fields being $L_{y}$ periodic, all non-zonal fluctuations are also $L_{x}$ periodic. However, this is not the case, for example, for the zonal density gradient $\partial _{x}\overline {n}$ , which becomes discontinuous at the end points (orange dots), even though the zonal density perturbation can be made periodic, as can be seen in figure 1. However, to use the pseudo-spectral algorithm, we need to compute its fast Fourier transform, since it appears in the nonlinear term $\partial _{x}\overline {n}\partial _{y}\phi$ .

To use fast Fourier transforms while maintaining non-periodic boundary conditions, we need to design a buffer zone at both ends of the $x$ axis to enforce the desired radial periodicity. For this purpose, we choose to restrict the radial domain $[0,L_{x}]$ to a smaller range $[x_{b1},x_{b2}]$ , with $0\lt x_{b1}\lt x_{b2}\lt L_{x}$ (see figure 2), which describes the physical region in which the solution is considered valid, and construct buffer zones, in which the profile is modified slightly to get the continuity of the first few derivatives of the zonal density profile (the continuity of the second derivative is needed to compute the neoclassical diffusion term).

Figure 2. Decomposition of the radial profile $n_{r}$ (blue) into its linear component $n_{lin}$ (dashed green) and its zonal fluctuating component $\overline {n}$ (dash-dotted red). The background density gradient is now computed between $x_{b1}$ and $x_{b2}$ , and the buffer zones are indicated by the dashed lines.

In the buffer region, the system has no ‘physical meaning’ and the profiles are not realistic, as a result of imposing the perturbation to match between the two ends of the domain. To avoid that these unrealistic profiles generate spurious transport, we also need to strongly damp fluctuations in those regions, which is achieved by the penalisation method.

When decomposing the radial profile, instead of computing the background density gradient using the endpoints of the domain $x=0$ and $x=L_{x}$ , we now take it between the physical end points $x_{b1}$ and $x_{b2}$ :

(3.1) \begin{equation} \kappa (t)\equiv -\frac {n_{r}(x_{b2},t)-n_{r}(x_{b1},t)}{x_{b2}-x_{b1}}\,. \end{equation}

An example of this decomposition is shown in figure 2.

Note that $x_{b2}-x_{b1}$ is equal to the physical box size $L_{px}$ , while $L_{x}$ , which also includes the two buffer regions, actually represents the computational domain size.

3.2. Flattening the zonal density in the buffer

The zonal profile inside the buffer zone needs to be such that at least its first two derivatives are continuous at the end points (i.e. $0$ and $L_{x}$ , which is the same point in a periodic box), while keeping the zonal perturbation periodic. This can be achieved by interpolating the zonal density perturbation in the buffer so that it basically becomes flat. However, the interpolation is rather costly from a computational point of view and, therefore, a better solution is to multiply the profile by a ‘gate function’, which is exactly 1 outside the buffer zones and 0 at the very edges. To ensure the continuity of the higher derivatives at the ends, the function that is used should be quite flat when converging to 0 (and $L_{x}$ on the other side) and it should be at least twice differentiable, to allow the computation of derivatives of the zonal profiles.

To guarantee that, we can use a smooth transition function (Tu Reference Tu2011, Ch. 13) which can be constructed using the following function, defined for $z\in [0,1]$ as

(3.2) \begin{equation} g(z)\equiv \begin{cases} \exp \left (-1/z\right ), & z\gt 0,\\ 0, & z=0. \end{cases} \end{equation}

This function is infinitely continuous and differentiable in its domain of validity and all its derivatives are equal to zero at $z=0$ , which makes it rather flat in its vicinity. Using this function as the basis, we can construct a smooth transition function, defined on $[0,1]$ as

(3.3) \begin{equation} h(z)\equiv \frac {g(z)}{g(z)+g(1-z)}, \end{equation}

which goes from 0 at $z=0$ to 1 at $z=1$ , and has all its derivatives equal to zero at both ends. The function $h$ is illustrated in figure 3(a).

Figure 3. (a) Smooth transition function $h$ defined according to (3.3) on $[0,1]$ . (b) Smooth gate $\varPsi$ defined from (3.4) on $[0,L_{x}]$ , where the transition parts are shown in orange.

Scaling and translating this function to match the buffer region, we can construct the following ‘smooth gate’ function $\varPsi$ :

(3.4) \begin{equation} \varPsi (x)\equiv \begin{cases} 0, & 0\leqslant x\leqslant x_{m_{\ell }},\\ h\left (\dfrac {x-x_{m\ell }}{x_{m1}-x_{m\ell }}\right ), & x_{m\ell }\lt x\lt x_{m1},\\ 1, & x_{m1}\leqslant x\leqslant x_{m2},\\ h\left (\dfrac {x_{mr}-x}{x_{mr}-x_{m2}}\right ), & x_{m2}\lt x\lt x_{mr},\\ 0, & x_{mr}\leqslant x\leqslant L_{x}, \end{cases} \end{equation}

which is exactly 0 in $[0,x_{m\ell }]$ and $[x_{mr},L_{x}]$ (close to the edges), exactly 1 on $[x_{m1},x_{m2}]$ , and which transitions smoothly for all its derivatives between those intervals. The gate function is shown in figure 3(b).

Notice that the intervals around the endpoints where the gate function is exactly zero are introduced for numerical reasons: since the $x$ axis is discretised in numerical simulations, it is better to ensure that the derivatives are exactly zero on a small portion of the $x$ grid, rather than just at the very ends. These regions $[0,x_{m\ell }]$ and $[x_{mr},L_{x}]$ are exagerated in figures 3 and 4 for visual purposes. In practice, we will choose $x_{m\ell }$ and $x_{mr}$ such that they correspond to some of the very first and last points of the grid (obviously within the buffer region), respectively (see § 4.1 for an actual application). In the following, we actually use the radial extension of the transition

(3.5) \begin{equation} \delta x_{m}\equiv x_{m1}-x_{m\ell }=x_{mr}-x_{m2} \end{equation}

as a parameter, and infer the individual values of $x_{m\ell }$ and $x_{mr}$ from that.

Figure 4. Top plot: modified zonal density perturbation $\overline {n}_{matched}$ (red), which is now periodic and flat at both ends. The corresponding matched radial profile $n_{r,matched}$ is shown in blue. The orignal radial profile $n_{r}$ (dashed light-blue), linear profile $n_{lin}$ (green) and zonal profile $\overline {n}$ (dashed dark-red) from figure 2 are also shown. Bottom plot: smooth-gate function $\varPsi$ .

Now, we can use the gate function $\varPsi$ to modify the zonal density profile inside the buffer zones, and match its values and derivatives at both edges. The simplest solution would be to just multiply the zonal density perturbation by $\varPsi$ to ensure the periodicity of $\overline {n}$ and all its derivatives. However, the values of the density perturbation, as computed from the decomposition, can be very different from zero at $x=0$ and $x=L_{x}$ , and forcing them to be zero could introduce strong gradients in the buffer zones which might trigger numerical instabilities (even though they would be suppressed in the buffer). To mitigate that effect, we first shift the zonal perturbation by an offset:

(3.6) \begin{equation} n_{off}\equiv \dfrac {1}{4}\left [\overline {n}(0)+\overline {n}(x_{m1})+\overline {n}(x_{m2})+\overline {n}(L_{x})\right ], \end{equation}

corresponding to the mean value between both ends of the profile and its value at $x=x_{m1}$ and $x=x_{m2}$ . Then, we apply the gate function and shift back the profile by $n_{off}$ :

(3.7) \begin{equation} \overline {n}_{matched}\equiv \varPsi (x)\left (\overline {n}-n_{off}\right )+n_{off}\,. \end{equation}

Note that this does not change the profile within the interval $[x_{m1},x_{m2}]$ . The matched zonal perturbation profile is shown on figure 4 (red line). We also show the corresponding ‘matched’ radial profile (blue line), which is the combination of the matched zonal perturbation and the linear profile,

(3.8) \begin{equation} n_{r,matched}\equiv n_{lin}+\overline {n}_{matched}, \end{equation}

actually ‘seen’ by the system.

Notice that we initially took $(x_{m1},x_{m2})=(x_{b1},x_{b2})$ , but this introduced strong density gradients very close to the physical space $[x_{b1},x_{b2}]$ , where the fluctuations are not completely zero. Instead, taking $x_{m1}$ and $x_{m2}$ well inside the buffer zone allows us to avoid such effects.

In practice, we subtract the mean value $n_{m}=\langle \overline {n}_{matched}\rangle _{x}$ from the matched zonal profile, as already discussed in § 2.2 (see (2.9)).

3.3. Penalisation method

Having constructed a profile whose zonal perturbation is periodic and whose derivative falls to zero at the endpoints, we now want to impose some constraints on both the radial profiles and turbulent fluctuations inside the buffer zone, to enforce the desired boundary conditions and avoid spurious fluctuations inside the edge regions. To this end, we use the volume penalisation approach, introduced by Angot et al. (Reference Angot, Bruneau and Fabrie1999) for Navier–Stokes simulations with a bounded, non-periodic domain, and used later for pseudo-spectral simulations of hydrodynamic (Schneider & Farge Reference Schneider and Farge2005) and magnetohydrodynamic (Schneider et al. Reference Schneider, Neffaa and Bos2011) turbulence.

The general idea is to impose a strong friction term

(3.9) \begin{equation} -\mu H(x)(A-A_{buff}) \end{equation}

to force the field $A$ to rapidly decay to a given boundary function $A_{buff}$ inside the buffer zone, represented by the mask function $H$ , which is 0 in the physical domain and 1 inside the buffer. To underline that it is large, the friction coefficient $\mu$ is sometimes represented as $1/\varepsilon$ , where $\varepsilon$ stands for some permeability coefficient, for the case when $A$ is the velocity field. This corresponds to building a thin artificial boundary layer around the buffer zone, in which the fields are forced to rapidly converge towards their edge value.

Using (3.9), and (2.11) and (2.15a ), we can write the penalised form of the equations that we solve numerically:

(3.10a) \begin{align} \partial _{t}\varOmega +[\phi ,\varOmega ] & =C(\widetilde {\phi }-\widetilde {n})+\nu {\nabla} ^{2}\widetilde {\varOmega }-\mu \boldsymbol{\nabla }\times \left (H(x)\boldsymbol{v}_{E\times B}\right ),\\[-10pt]\nonumber \end{align}

(3.10b) \begin{align} \partial _{t}N+[\phi ,N] & =C(\widetilde {\phi }-\widetilde {n})+D{\nabla} ^{2}\widetilde {n}+S_{n}(x,t)-\mu H(x)\left (N-n_{r,buff}\right ), \end{align}

where the last two terms in both equations correspond to the penalisation. Notice that in (3.10a ), we impose the fluid velocity $\boldsymbol{v}_{E\times B}=-\partial _{y}\phi \boldsymbol{e}_{x}+\partial _{x}\phi \boldsymbol{e}_{y}$ to be zero inside the buffer zone, corresponding to no-slip boundary conditions, and since the equation is actually for the vorticity $\varOmega =\boldsymbol{\nabla }\times \boldsymbol{v}_{E\times B}$ , it is the curl of that penalisation condition (projected onto $\boldsymbol{e}_{z}$ ) that appears in this equation (Schneider & Farge Reference Schneider and Farge2005; Schneider et al. Reference Schneider, Neffaa and Bos2011). For the density (3.10b ), the total density $N$ is forced to be equal to a given boundary function $n_{r,buff}(x,t)$ in the buffer, which is only a function of $x$ and time, so that the turbulent, non-zonal fluctuations are forced to be zero in the buffer. The choice for $n_{r,buff}$ will be discussed further in § 3.4 about boundary conditions on the radial profile.

We can thus write the penalised equations for the turbulent fluctuations (using (2.15a ) and (2.15b )),

(3.11a) \begin{align} \partial _{t}\widetilde {\varOmega } & +[\phi ,\widetilde {\varOmega }]-\langle [\phi ,\widetilde {\varOmega }]\rangle _{y}=C(\widetilde {\phi }-\widetilde {n})+\nu \boldsymbol{\nabla} ^{2}\widetilde {\varOmega }-\mu \boldsymbol{\nabla }\times \left (H(x)\hat {z}\times \boldsymbol{\nabla }\widetilde {\phi }\right ), \end{align}

(3.11b) \begin{align} \partial _{t}\widetilde {n} & +\left (\kappa (t)-\partial _{x}\overline {n}\right )\partial _{y}\widetilde {\phi }+\partial _{x}\overline {\phi }\partial _{y}\widetilde {n}+\widetilde {[\phi ,n]}=C(\widetilde {\phi }-\widetilde {n})+D\boldsymbol{\nabla} ^{2}\widetilde {n}-\mu H(x)\widetilde {n}, \end{align}

and for the radial profile of the velocity $u_{r}$ , which is just the zonal velocity $\overline {v}_{y}$ here,Footnote ² and the radial density $n_{r}$

(3.12a) \begin{align} \partial _{t}\overline {v}_{y} & =\langle \widetilde {\varOmega }\partial _{y}\widetilde {\phi }\rangle _{y}-\mu H(x)\overline {v}_{y}\\[-10pt]\nonumber \end{align}

(3.12b) \begin{align} \partial _{t}n_{r} & =\partial _{x}\langle \widetilde {n}\partial _{y}\widetilde {\phi }\rangle _{y}+S_{n}(x,t)+D_{0}\partial _{x}^{2}n_{r}-\mu H(x)\left (n_{r}-n_{r,buff}\right )\,. \end{align}

The mask function $H(x)$ , which should quickly (but smoothly) transition from 0 to 1 as we enter the buffer zone, is taken as the smooth gate function $\varPsi$ defined previously in (3.4). More precisely, we will use the following mask function:

(3.13) \begin{equation} H(x)=1-\varPsi _{mask}(x)\,. \end{equation}

Note that $\varPsi _{mask}$ is parametrised differently from the function used to make zonal profiles periodic in § 3.2. As illustrated in figure 5, the parameters are chosen so that fluctuations are zero before the radial profile has any modification. In other words the fluctuations do not ‘see’ the modified parts of the radial profiles.

Figure 5. Combined plot of $\psi _{mask}(x)=1-H(x)$ (blue) and the smooth-gate function $\varPsi _{smooth}(x)$ (dashed red). Their radial extensions $\delta x_{b}$ and $\delta x_{m}$ are also shown (respectively blue and red arrows).

3.4. Boundary conditions

3.4.1. Boundary profile in the buffer zone

First, we specify the boundary radial profile $n_{r,buff}(x,t)$ that we impose in the buffer zone. Since the fluctuations are not supposed to penetrate that region, we argue that the density profile should keep its initial shape $n_{r}(x,t=0)$ . However, since turbulence creates a net outward radial particle flux, the total number of particles can vary in time in the physical domain. This would create large gradients in the buffer zone if we kept the density in the buffer region really constant. To avoid this, we need the boundary profile to decrease or increase together with the value of $n_{r}(x,t)$ at $x_{b1}$ and $x_{b2}$ , where the buffers meet the physical domain. To achieve this, we propose the following boundary profile:

(3.14) \begin{equation} n_{r,buff}(x,t)=\begin{cases} n_{r}(x,t=0)-n_{r}(x_{b1},t=0)+n_{r}(x_{b1},t), & x\leqslant x_{b1},\\ n_{r}(x,t=0)-n_{r}(x_{b2},t=0)+n_{r}(x_{b2},t), & x\geqslant x_{b2}, \end{cases} \end{equation}

which will cause the boundary profiles to move together with $n_{r}(x_{b1},t)$ and $n_{r}(x_{b2},t)$ at these points, without changing their initial shape.

3.4.2. Boundary condition at $(x_{b1},x_{b2})$

As a first application of the code with non-trivial boundary conditions, we choose to fix the density profile at the right buffer, $x=x_{b2}$ (Dirichlet boundary condition), while letting it evolve freely at the left buffer, $x=x_{b1}$ , as a response to the turbulent particle flux (which is a loosely Neumann boundary condition). Note that we make this somewhat arbitrary choice because we are interested in studying the profile relaxation, or its evolution with a particle source localised at the inner boundary, which is more relevant if the profile is fixed at the outer boundary. It also allows us to demonstrate the two kinds of boundary conditions that can be used.

In particular, imposing the Dirichlet boundary condition $\partial _{t}n_{r}(x_{b2},t)=0$ is not straightforward: even though we impose a boundary profile for $n_{r}$ inside the buffer, the condition is only really met where the mask function $H(x)$ is close to $1$ . Therefore, in the vicinity of $x_{b2}$ in the buffer, where $H(x)\ll 1$ (see figure 5), the radial profile can still evolve freely because the effect of penalisation is quite weak in this region. To impose the condition exactly at $x_{b2}$ , but in a smooth way around this point, we apply an artificial Gaussian sink centred at the outer boundary point. Indeed, we cannot impose a boundary condition only at $x_{b2}$ , because such a singularity would trigger numerical errors when computing fast Fourier transforms.

This artificial Gaussian sink is taken as

(3.15) \begin{equation} S_{b2}(x,t)=-\partial _{t}n_{r}(x_{b2},t)\exp \left [-(x-x_{b2})^{2}/(2\sigma _{S_{b}}^{2})\right ], \end{equation}

where $\partial _{t}n_{r}(x_{b2},t)$ is computed from (3.12b ). The sink is centred on $x_{b2}$ and imposes $\partial _{t}n_{r}(x_{b2},t)=0$ . Since it has some finite (small) width $\sigma _{S_{b}}$ , its effect extends into the physical space around $x_{b2}$ , which means that it modifies the particle budget inside the physical domain and acts as a volumetric particle sink. The particle budget becomes

(3.16) \begin{equation} d_{t}N_{\varphi }=\int _{x_{b1}}^{x_{b2}}\partial _{t}n_{r}\,\text{d}x=\varGamma _{n}(x_{b1},t)-\varGamma _{n}(x_{b2},t)+P_{b2}(t), \end{equation}

with $N_{\varphi }$ the total density inside the physical domain $[x_{b1},x_{b2}]$ , $\varGamma _{n}$ the total particle flux (turbulent and potentially diffusive if the neoclassical diffusion $D_{0}$ is also included) and $P_{b2}=\int _{x_{b1}}^{x_{b2}}S_{b2}(x,t)\,\text{d}x$ is the integral over the physical domain of the sink $S_{b2}$ localised at $x_{b2}$ .

Note that this is a particular choice for the present study and that the boundary conditions can easily be modified for other applications of the code.

3.5. Limits of the penalisation method

First of all, notice that, by construction, turbulent (zonal and non-zonal) perturbations are taken to be zero mean-valued. However, since we impose the radial profile of the poloidal velocity $u_{r}$ to be zero in the buffer zone, it may acquire a non-zero mean value resulting in $\overline {v}_{y}=u_{r}-\langle u_{r}\rangle _{x}\neq u_{r}$ . We account for this small non-zero mean value in the numerical implementation, which can be seen as a Doppler shift (and which remains several orders of magnitude lower than the typical amplitude of $u_{r}$ in the saturated state). This makes it so that the zonal velocity perturbation is not exactly zero in the buffer zone, but it is actually equal to $-\langle u_{r}\rangle _{x}$ . This is not really a problem, but a feature of the penalisation with the particular kind of boundary condition that is used.

Second, the actual effectiveness of the penalisation is related to the value of the penalisation coefficient $\mu$ in such a way that the higher the coefficient $\mu$ , the closer we are to a system with actual impermeable boundaries. However, since $\mu \lt \infty$ , there is still some permeability which means that some fluctuations manage to penetrate into the buffer region. Nevertheless, we observe in simulations that the turbulence level in the buffer zone is still 2–3 orders of magnitude below that of what one finds in the physical domain and no coherent information seems to pass from one buffer to the other. This can, in principle, be improved by increasing the numerical value of the penalisation coefficient, but experience shows that the improvement remains limited, and it causes the adaptive time stepping algorithm that we use to slow down considerably.

The mask function can also create a second instability induced by the penalisation term $-\mu \boldsymbol{\nabla }\times (H(x)\hat {z}\times \boldsymbol{\nabla }\widetilde {\phi })$ from (3.11a ). Indeed this term can be decomposed into a penalisation term on the vorticity $-\mu H(x){\nabla} _{\perp }^{2}\widetilde {\phi }$ and a second term $-\mu \boldsymbol{\nabla }H\times (\hat {z}\times \boldsymbol{\nabla }\widetilde {\phi })$ , which features the gradient of the mask function and which will appear as an additional term in the dissipative drift-wave dispersion relation when kept, in contrast to Navier–Stokes, where there is no linear instability. Note that this is not a numerical artefact per se: this would probably happen if we had actual walls around the system as implied by the penalisation scheme. However, it is unwelcome here because the buffer zones are not supposed to represent actual walls and it shows that using the penalisation method may cause the dynamics close to the boundaries to be affected in specific ways, which may be different from what is intended. Fortunately, this instability remains of very low amplitude, confined to the vicinity of the buffer zone, and is completely suppressed once ‘real’ turbulence reaches the buffer and covers it (see the discussion in § 4.3.1 on an actual DNS). Again, this effect can be reduced by making the buffer regions larger, so that the local gradients of the mask function become smaller. However, this is a trade-off, since in doing so, one looses a bigger part of the computational domain to buffers.

4.1. Numerical set-up

To further study the dynamics of the flux-driven Hasegawa–Wakatani system accurately, we perform high-resolution ( $4096\times 4096$ ) numerical simulations, for which the physics and penalisation parameters are given in table 1. The parameters that determine the buffer and smoothing functions correspond to integer multiples of the unpadded grid resolution $\Delta x=L_{x}/N_{x}$ , where $N_{x}=2\lfloor 4096/3\rfloor =2730$ . The initial value for the profile is given by (4.1) as before.

Table 1. Physical and penalisation parameters for the $4096\times 4096$ padded simulation. The initial profile is taken as (4.1).

Note that for this high-resolution simulation, the buffer zone corresponds to less than $7\,\%$ of the total computational domain size. In contrast, in the $1024\times 1024$ simulation shown in figure 6, the buffer zone covers $25\,\%$ of the computational domain size (using the same integer multiples for setting the buffer zone).

The time evolution of energy, flux and density at the left boundary $n_{r}(x_{b1},t)$ are shown in figure 18 in Appendix B.

4.2. Early times: pair production

First, we detail the flattening of the initial profile as a result of positive blobs of density $\delta n_{+}$ travelling outwards and negative holes $\delta n_{-}$ travelling inwards. This can be seen as a consequence of the modulation of the linear drift-wave instability, which occurs mostly where the profile is the steepest (highest local $\kappa$ ). Although the waveform of this instability is mainly poloidal (i.e. $k_{y}=0$ modes have zero growth rate), the nonlinear interaction of two linearly unstable modes can form a structure with a finite radial scale through a mechanism akin to modulational instability.

The evolution of the profile is shown in figure 7(a) at different time steps (the initial profile is the black dashed line). The flattening is centred at $x=x_{a}\approx 53$ , around which the modulation grows, making the profile more steep at the edges. The perturbation with respect to the initial profile

(4.2) \begin{equation} \delta n_{r}\equiv n_{r}-n_{r}(t=0) \end{equation}

is shown in figure 7(b), where we clearly see the formation of a negative ‘dip’ on the left, corresponding to the holes $\delta n_{-}$ , and a positive bump on the right, corresponding to the positive blobs $\delta n_{+}$ , both moving away from $x=x_{a}$ and, as a result, flattening the profile. The growth of this flattening perturbation, which we may call an avalanche from the point of view of the profile, is initially smooth, but becomes chaotic towards the end $t=30.6$ (violet line). At this point, the system really enters the turbulent regime, as the blobs become unstable themselves, forming smaller scale structures, as can be seen in figures 6(c) and 6(g).

Figure 7. Flattening of the initial density profile in the early times. (a) Density profile restrained on the interval $x\in [30,75]$ at different time steps. The initial profile is in black dashed lines. (b) Perturbation $\delta n_{r}=n_{r}-n_{r}(t=0)$ from the intial profile.

To quantify the dynamics, we can look at the time evolution of the root mean square (r.m.s.) of the deviation from the initial profile

(4.3) \begin{equation} \langle \delta n_{r}\rangle _{rms}\equiv \left (\frac {1}{x_{2}-x_{1}}\int _{x_{1}}^{x_{2}}\delta n_{r}^{2}\,{\rm d}x\right )^{1/2}, \end{equation}

where we take $x_{1}=30$ and $x_{2}=75$ since it is sufficient to look only in the vicinity of the avalanche. We can also look at the position $X_{\pm }$ of the positive and negative perturbations $\delta n_{\pm }$ that we observe in figure 7(b). We define those as the barycentre of the squared perturbations:

(4.4) \begin{equation} X_{\pm }(t)\equiv \frac {\int _{x_{1}}^{x_{2}}x\delta n_{r}^{2}\varTheta \left (\pm \delta n_{r}\right )\,{\rm d}x}{\int _{x_{1}}^{x_{2}}\delta n_{r}^{2}\varTheta \left (\pm \delta n_{r}\right )\,{\rm d}x}, \end{equation}

where $\varTheta$ is the Heaviside step function. This definition is similar to tracking the position of the maxima, but it has the advantage of behaving more smoothly.

In figure 8, we show the r.m.s. value $\langle \delta n_{r}\rangle _{rms}$ of the perturbation (panel a), and the position $X_{\pm }$ of the bump and the hole (panel b) as functions of time. The time evolution of the perturbation can be decomposed in two phases. First, there is a linear phase, where the perturbation grows exponentially and the bump and the hole stay almost at the same distance from each other, although we can notice that they slowly drift towards the left (in other simulations, they may not move at all). Then, for $t\geqslant 27$ , the amplitude of the perturbation saturates, which is associated with the positive perturbation moving radially outwards while the negative perturbation moves inwards, at roughly constant speeds. Finally, towards the end, the perturbation starts to be chaotic and the blobs become unstable themselves.

Figure 8. (a) Time evolution of the r.m.s. value $\langle \delta n_{r}\rangle$ of the density profile perturbation. The slope of the growth of the perturbation in the linear phase is compared with the sum of the growth rates of the most unstable mode and side-band $\gamma _{k}+\gamma _{p}$ . (b) Positions $X_{\pm }$ of the density perturbations $\delta n_{\pm }$ as functions of time, shown respectively in red and blue dashed lines. The contour plots of the density perturbation $\delta n_{r}$ are also shown. The perturbation wavelength $\lambda _{q}=2(X_{+}-X_{-})\approx 23.0$ is estimated at $t\approx 20$ (orange arrow).

To explain the formation of the perturbation, which is a radial structure, we can assume that it results from triadic interactions between the most unstable mode $k=(0,k_{y0})$ , a radial mode $q=(q,0)$ , with $q$ corresponding to the wavenumber of the perturbation, and two side-bands $p_{\pm }=(\pm q,k_{Y0})$ , which statisfy the triadic interaction $\boldsymbol{k}-\boldsymbol{p}_{\pm }\pm \boldsymbol{q}=\boldsymbol{0}$ . Taking the Fourier transform of (2.12), and restraining the nonlinear term to only those, yields

(4.5) \begin{equation} \partial _{t}\overline {n}_{q}=qk\left (\phi _{k}^{*}n_{p_{+}}-\phi _{p_{+}}n_{k}^{*}-\phi _{k}n_{p_{-}}^{*}+\phi _{p_{-}}^{*}n_{k}\right )\,. \end{equation}

Assuming $\phi _{k,p_{\pm }},n_{k,p_{\pm }}\propto \exp [-i\omega _{k,p_{\pm }}t]$ , where $\omega _{k}=\omega _{k,r}+i\gamma _{k}$ is the eigenvalue associated with the unstable mode ( $\gamma _{k}\gt 0$ ) of the drift-wave instability, we can write

(4.6) \begin{equation} \overline {n}_{q}=\left (A_{k,q}^{+}e^{-i\delta \omega t}+A_{k,q}^{-}e^{i\delta \omega t}\right )e^{(\gamma _{k}+\gamma _{p})t}, \end{equation}

where $\delta \omega =\omega _{p,r}-\omega _{k,r}$ and $A_{k,p}^{\pm }$ are complex coefficients depending on the initial phase and amplitudes of the modes $k$ and $p_{\pm }$ . This suggests that the pertubation is expected to grow at a rate $\gamma _{q}^{\delta n}=\gamma _{k}+\gamma _{p}$ initially.

Note that from figure 8(b), we can estimate the perturbation wavelength to be $\lambda _{q}=2(X_{+}-X_{-})\approx 23.0$ (at $t\approx 20$ , when the distance between the bump and the hole is still constant), which yields a wavenumber $q\approx 0.27$ . The most unstable wavenumber can be calculated using the dispersion relation (see Appendix C and Gürcan et al. (Reference Gürcan, Anderson, Moradi, Biancalani and Morel2022) for the analytical expression and discussions) with $C=0.05$ and $\kappa _{\ell }=5.0$ , and yields $k_{y0}\approx 0.39$ . We can then compute the perturbation growth rate $\gamma _{k}+\gamma _{p}\approx 0.73$ , which corresponds to the dashed black line in figure 8(b), and which shows a very good agreement between this value and the actual growth of $\delta n_{r}$ in the simulation.

The mechanism that determines the wavenumber $\lambda =2\pi /q$ of the perturbation is yet to be found, but it might be similar to the one that determines the typical size of the zonal flows, which also remains unknown in the case of the Hasegawa–Wakatani system. Notice also that the pair formation is only observed if the sharp gradient region is sufficiently extended compared with $\lambda =2\pi /q$ . If it is too narrow, we do not see the formation of these large-scale corrugations, but rather smaller, less coherent structures.

4.3. Relaxation and spreading

4.3.1. Observations of profile relaxation and turbulent spreading in the DNS

To understand how the system continues to evolve, we consider the spreading of the turbulent region together with the relaxation of the radial density profile over longer time scales. In figure 9(a), we show the density profile $n_{r}$ at different time steps, starting from $t=30$ (dark blue), when the initial pair consisting of a bump and a hole collapses (the initial density profile is shown in black dashed lines), until $t=1000$ (brown) when turbulence reaches the right buffer and zonal flows start to accumulate and dominate the system. A time animation of both density and zonal velocity radial profiles is available on the repository page (Guillon Reference Guillon2025a , Movie 2). It can be seen that the initial sharp gradient gradually relaxes in time, as it extends to the right into the initially flat part of the profile. Spreading to the right slows down and stops when it reaches the right buffer, at which point the system is completely covered by turbulence, as can be seen in figures 6(d) and 6(h). Notice that the profile completely looses its original shape and is very close to a straight line as a result of the turbulent particle flux, flattening all the steep, large amplitude initial perturbations.

Figure 9. (a) Radial density profile $n_{r}$ at different time steps (the dashed line correspond to the initial profile). (b) Radial turbulent kinetic energy profile $\overline {\mathcal{K}}$ at the same time steps, in $y$ semilog plot. Circles correspond to the right-most position at which the radial kinetic energy $\overline {\mathcal{K}}$ is equal to $1\,\%$ of its maximum at a given time step (4.9), which we use as a proxy for the turbulent front position $X_{f}(t)$ .

We can also study the radial profile of the turbulent intensity and its evolution in time. To do so, we define the radial turbulent kinetic energy as

(4.7) \begin{equation} \overline {\mathcal{K}}(x)\equiv \langle \widetilde {v}_{x}^{2}+\left (\widetilde {v}_{y}+\overline {v}_{y}\right ){}^{2}\rangle _{y}, \end{equation}

which is the poloidal average of the perturbed kinetic energy $e_{\mathcal{K}}(x,y)=\widetilde {v}_{x}^{2}+(\widetilde {v}_{y}+\overline {v}_{y}){}^{2}$ , where we include both zonal and non-zonal perturbations. The time evolution of the radial profile of kinetic energy is shown in figure 9(b), at the same time steps with the density profile in panel (a). The first time step $t=30$ (dark blue) shows that the turbulent energy is peaked around $x=x_{a}$ , where the profile is the steepest. In fact, we could obtain this shape by solving the equation

(4.8) \begin{equation} \partial _{t}\overline {\mathcal{K}}=2\gamma _{max}(\partial _{x}n_{r})\overline {\mathcal{K}}, \end{equation}

where $\gamma _{max}(\partial _{x}n_{r})$ is the maximum linear growth rate computed by taking $\kappa (x)=-\partial _{x}n_{r}$ at each point $x$ and is thus a function of the radial position. This is discussed in more detail in § 4.3.2, where we compare the DNS results with those of a 1-D reduced model of transport/turbulence spreading. After the initial growth, the kinetic energy peak saturates and spreads in the radial direction, together with the relaxing radial density profile, decreasing in amplitude as the local free energy source is lowered because the local profile gradient decreases. At $t=1000$ (brown), turbulence completely covers the radial domain, and we start to see some large-scale peaks at $x\approx 40$ and $x\approx 185$ , which corresponds to the formation of zonal flows. Notice that for some intermediary time steps (e.g. $t=100$ , in light-blue), we can observe traces of the numerical instability due to the mask function gradient in the vicinity of the right buffer, as discussed in § 3.5. As mentioned, this instability is 2–3 orders of magnitude lower than the turbulence level and becomes completely irrelevant once the turbulence reaches the right buffer.

To investigate how the density profile relaxation/turbulence spreading front moves, we need a proxy for the position of the front so that we can study its propagation. Ideally, the position of this front should correspond to the point at which the density profile is turbulent on the left and still at rest on the right. To this end, we use as a proxy the last radial point where the turbulent kinetic energy exceeds a certain threshold, which we take as $1\,\%$ of the maximum energy at each time step:

(4.9) \begin{equation} X_{f}(t)\equiv \max _{x}\big\{ x\;:\;\overline {\mathcal{K}}(x,t)\geqslant 0.01\times \max _{x}\big [\overline {\mathcal{K}}(x,t)\big ]\big\} \,. \end{equation}

The position of $X_{f}$ at various time steps is shown in figure 9 (coloured circles). We see that the proxy provides a good estimate of the position of the front.

In figure 10(a), we show the time evolution of the turbulent front position $X_{f}(t)$ (dashed red line), on top of the spatiotemporal evolution of the radial profile of kinetic energy $\overline {\mathcal{K}}(x,t)$ . We find again that the localised growth of the turbulent energy around $x=x_{a}$ spreads quickly until $t=100$ and then slows down but continues spreading at a slower rate, while decreasing in amplitude. Note that we do not plot $X_{f}$ for $t\lt 30$ because the definition of the proxy (4.9) fails when the system is still in the linear instability phase, concentrated at $x_{a}$ . At $t\approx 600$ , we start to see the formation of zonal flows, evidenced by the yellow horizontal stripes in the spatiotemporal evolution of $\overline {\mathcal{K}}$ .

Figure 10. (a) Time evolution of the turbulent front position $X_{f}(t)$ (dashed red line), on top of the spatiotemporal evolution of the radial kinetic energy profile $\overline {\mathcal{K}}(x,t)$ . (b) Time evolution of the left mean density gradient, defined by (4.11). (c) Time evolution of the corresponding maximum growth rate, defined by (4.11). Both mean gradient and growth rate are fit with power laws of the time (red and green dashed lines). The vertical grey line denotes the time $t_{1}\approx 800$ at which the front reaches the right buffer.

Using the turbulent front position, we can track the time evolution of the mean density gradient to its left, $\kappa _{\ell }$ , which corresponds to the slope of the high gradient part of the profile, by using

(4.10) \begin{equation} \kappa _{\ell }(t)\equiv \frac {n_{r}(x_{b1},t)-n_{r}(X_{f}(t),t)}{X_{f}(t)-x_{b1}}, \end{equation}

which is defined using the difference between the left buffer position $x_{b1}$ and the right turbulent front position $X_{f}(t)$ . Note that maybe we should rather define a left turbulent front position and use it instead of $x_{b1}$ . However, since the left front reaches the left buffer very quickly, it is practically equivalent to using $x_{b1}$ . The time evolution of $\kappa _{\ell }(t)$ is shown in figure 10(b). After the initial collapse, the mean left gradient decreases according to a time power law $t^{-0.46}$ (red dashed line) while turbulence spreads to the right. When the turbulent front reaches the right buffer ( $t\geqslant 800$ ), $\kappa _{\ell }$ becomes the global density gradient of the system and continues to decrease according to a slower power law $t^{-0.13}$ (green dashed line), due to the difference of particle flux between both ends of the physical domain, as there is no particle source in these simulations to compensate for the flux difference.

We can also compute the maximum growth rate in the left part associated with $\kappa _{\ell }(t)$ ,

(4.11) \begin{equation} \gamma _{\ell }(t)\equiv \gamma _{max}\left (\kappa _{\ell }(t)\right ), \end{equation}

which corresponds to the energy injection rate in the turbulent region. The time evolution of $\gamma _{\ell }(t)$ , shown in figure 10(c), is very similar to that of $\kappa _{\ell }$ , and also follows two different time power laws before and after the front reaches the buffer, though the power exponents slightly differ from that of $\kappa _{\ell }(t)$ . From the figure, we see that the energy injection rate decreases with the relaxation of the profile. This can actually be used as an estimate of the slow turbulent spreading velocity, which is rather low for this system, as discussed in the next section.

4.3.2. One-dimensional reduced model for coupled profile relaxation/spreading

To understand the key mechanisms that play a role in the nature of profile relaxation/turbulence spreading that we observe in the DNS without explicit source terms, we use a 1-D reduced model which features an equation on the radial turbulent kinetic energy $\overline {\mathcal{K}}$ coupled with a transport equation for the density profile $n_{r}$ . The first equation is

(4.12) \begin{equation} \partial _{t}\overline {\mathcal{K}}=2\gamma _{max}(x,\partial _{x}n_{r})\overline {\mathcal{K}}-\beta _{NL}\overline {\mathcal{K}}^{2}+\chi _{\mathcal{\mathcal{K}}}\partial _{x}\left (\overline {\mathcal{K}}\,\partial _{x}\overline {\mathcal{K}}\right ), \end{equation}

where on the right-hand side, the first term corresponds to the local linear maximum growth rate (i.e. the growth rate of the most unstable mode $k_{y0}$ ), computed using the radial density gradient $\partial _{x}n_{r}$ at each point as the ‘local’ background density gradient. The second term corresponds to a nonlinear saturation due to eddy viscosity, determined by $\beta _{NL}$ , so that we would have a steady-state $\overline {\mathcal{K}}(x)=2\gamma _{max}(x)/\beta _{NL}$ without turbulent spreading. The latter is in turn represented by the last term on the right-hand side, which is a nonlinear diffusion term, with a diffusion coefficient $\chi _{\mathcal{K}}\overline {\mathcal{K}}$ that is proportional to the turbulent kinetic energy itself. The second equation, for the radial density profile, can be written as

(4.13) \begin{equation} \partial _{t}n_{r}=D_{n}\partial _{x}\left (\overline {\mathcal{K}}\,\partial _{x}n_{r}\right ), \end{equation}

where the turbulent flux is also assumed to have the same nonlinear diffusivity as the turbulent kinetic energy, i.e. $\varGamma _{n}=-D_{n}\overline {\mathcal{K}}\,\partial _{x}n_{r}$ . The coefficients $\beta _{NL},\chi _{\mathcal{K}}$ and $D_{n}$ are the parameters of the model, which is very similar to previous reductions introduced to study profile relaxation and turbulence spreading together. Although here, (4.12) is on the radial turbulent kinetic energy $\overline {\mathcal{K}}$ , for some of these works, it is replaced by the evolution of the turbulent intensity $\sum _{k}|\phi _{k}|^{2}$ denoted by $I$ (Hahm et al. Reference Hahm, Diamond, Lin, Itoh and Itoh2004; Miki et al. Reference Miki, Diamond, Gürcan, Tynan, Estrada, Schmitz and Xu2012), $\varepsilon$ (Gürcan et al. Reference Gürcan, Diamond, Hahm and Lin2005) or $E$ (Garbet et al. Reference Garbet, Sarazin, Imbeaux, Ghendrih, Bourdelle, Gürcan and Diamond2007). Naulin et al. (Reference Naulin, Nielsen and Rasmussen2005) use an equation for a general definition of the turbulent energy, written as $\epsilon$ .

To compare the 1-D reduced model with DNS, we incorporate the boundary conditions of the penalisation method in (4.12) and (4.13) (see Appendix D). Note that some refined model (Miki et al. Reference Miki, Diamond, Gürcan, Tynan, Estrada, Schmitz and Xu2012) can include an additional saturation term due to zonal flow shearing, as well as an equation on the zonal energy. However, as will be discussed in § 4.4, the fraction of zonal energy remains rather low during the spreading phase, and the two-equations system (4.12) and (4.13) seemed to be sufficient to reproduce the turbulent front spreading, in our case. Additional refinements include a mean $E\times B$ poloidal flow shearing effect, along with neoclassical diffusion and flow terms, which are out of the scope of the present study.

In figure 11, we show the comparison of the turbulent front position $X_{f}(t)$ (panel a) and velocity $v_{X_{f}}(t)=d_{t}X_{f}$ (panel b) measured using both DNS and the 1-D reduced model. The numerical simulation of the 1-D reduced model is performed using $\beta _{NL}=0.014,\chi _{\mathcal{K}}=0.1$ and $D_{n}=0.1$ , which yields satisfactory agreement with the DNS, even though a formal optimisation study, which is left for future work, could be used to improve this agreement. The turbulent front position $X_{f}(t)$ is measured using the proxy (4.9), from which we subtracted its initial value at $t_{0}=30$ (before this time, the proxy measurement fails because we are in the linear phase). The velocity $v_{X_{f}}(t)=d_{t}X_{f}$ is measured in both cases by differentiating $X_{f}$ after applying a Savitzky–Golay filter (Savitzky & Golay Reference Savitzky and Golay1964) with a window of 1000 time points and a third-degree polynomial, to filter out the short fluctuations of the front position (which arise from the front position measurement using the proxy and from fast turbulent fluctuations), which correspond to abrupt and noisy spikes in the velocity signal. That way, we can focus on the mean velocity of the front which evolves on a slower time scale than that of turbulence.

Figure 11. Comparison of the turbulent spreading between the DNS (blue) and the 1-D reduced model (red). (a) Log–log plot of the turbulent front $X_{f}(t)$ computed using the proxy (4.9), from which the initial front position $X_{f}(t_{0})$ has been subtracted. The green dashed line corresponds to the 1-D simulation without spreading (i.e. $\chi _{\mathcal{K}}=0$ ). The black dashed line corresponds to the slope of the diffusive spreading $X(t)\propto t^{1/2}.$ (b) Velocity of the turbulent front $v_{X_{f}}(t)=d_{t}X_{f}$ computed by differentiating $X_{f}(t)$ after applying a Savitzky–Golay filter. The simulations are compared with the theoretical estimation of the front velocity $\sqrt {D_{n}\gamma _{\ell }(t)}$ (dashed line) where $\gamma _{\ell }$ is computed using (4.11) and $D_{n}=0.1$ is the diffusion coefficient of the 1-D reduced density transport (4.13).

The comparison of the turbulent front behaviour between the DNS and the 1-D reduced model yields remarkably good agreement, especially for the long time evolution of the front position $X_{f}(t)$ , where it seems to follow some kind of power law with $X_{f}(t)\propto t^{a}$ , where $a\approx 0.4$ . This suggests that turbulence spreading is slightly subdiffusive in this particular example, since $X_{f}(t)$ evolves slower than $t^{1/2}$ , as shown by the black dashed line in panel (a). Moreover, in different simulations of the 1-D reduced model, we oberve that the exponent $a$ of the power law scaling increases as a function of $\chi _{\mathcal{K}}$ and $D_{n}$ (taken equal).

To underline the importance of the turbulent spreading term $\chi _{\mathcal{\mathcal{K}}}\partial _{x}(\overline {\mathcal{K}}\,\partial _{x}\overline {\mathcal{K}})$ , we performed a simulation of the 1-D model taking $\chi _{\mathcal{K}}=0$ , while keeping the same values for $\beta _{NL}$ and $D_{n}$ . The front propagation of this simulation without the spreading term corresponds to the green dashed line in figure 11(a). We see that in this case, the front does not propagate as far (two times less) as when turbulent spreading is included, and it really becomes diffusive in this case (i.e. evolves as ${t}^{1/2}$ ). To observe a similar front propagation, we need to multiply the diffusion coefficient $D_{n}$ approximately by 10. In particular, we loose the initial, brutal collapse of the profile ( $t=50{-}100$ in figure 10 a). Hence, even though turbulence spreading makes the front propagation slightly subdiffusive, its diffusion rate becomes much (almost 10 times) higher than profile relaxation on its own, suggesting that it plays a role in the front propagation in this system.

The front velocity closely matches between the DNS and the reduced model, although it fluctuates much more in the former due to the complete turbulence description (with eddies intermittently surging into the unperturbed region on the right), while it propagates more smoothly in the 1-D model where the spreading comes from the nonlinear diffusive term. We compare the numerical results to the theoretical estimate of the front velocity for a diffusive process $v_{front}=\sqrt {D_{n}\gamma _{\ell }(t)}$ (Gürcan et al. Reference Gürcan, Diamond, Hahm and Lin2005), where $\gamma _{\ell }(t)$ is the maximum growth rate in the turbulent domain (measured in the DNS), and $D_{n}$ is the density coefficient from the reduced model. Although this estimation fails at the early times, where the velocity decreases very rapidly, it roughly matches both models at long time, as evidenced by the black line in figure 11(b).

The animation of the energy and density profiles for both DNS and the 1-D model is available at the repository page (Guillon Reference Guillon2025a , Movie 3), and shows that the profiles can be different to some extent between both simulations, even though the front evolution is very similar. This could be improved by optimising the parameters of the 1-D reduced model using the DNS to completely validate the comparison, along with adding the effect of turbulence damping by zonal flows shear, which is left for future work.

4.4. Freezing of the system by zonal flows

Last, we can discuss the final state of the system, after $t=t_{1}$ , at which the turbulent front reaches the right buffer. The density profile continues to relax until it stops and develops large-scale corrugations. These corrugations can be seen as rising trends and sharp drops in the zonal density profile $\overline {n}$ (blue line), shown in figure 12(a), together with smaller scale fluctuations. At the same time, we see the emergence of large-scale zonal flows (orange) which become stationary.

Figure 12. (a) Zonal density profile $\overline {n}$ (blue) and zonal velocity profile $\overline {v}_{y}$ (orange) at the last time step of the simulation. (b) Control parameter $C/\kappa$ computed using the left density gradient $\kappa _{\ell }(t)$ (blue) and the global mean gradient $\kappa$ (dashed green) as a function of time. The fraction of zonal kinetic energy $\varXi _{\mathcal{K}}$ is also shown in red. The vertical line marks the time $t_{1}\approx 800$ at which the turbulent front reaches the right buffer.

This is due to the fact that, since the turbulent front has reached the right buffer, the density gradient $\kappa _{\ell }$ , which now corresponds to the mean global background $\kappa$ of the system, is close to the critical value $\kappa _{c}=0.5$ so that $C/\kappa _{c}=0.1$ , at which the fixed gradient system normally transitions to a zonal flow dominated state (Numata et al. Reference Numata, Ball and Dewar2007; Grander et al. Reference Grander, Locker and Kendl2024; Guillon & Gürcan, Reference Guillon and Gürcan2025). This is illustrated in figure 12(b) where the ratios $C/\kappa _{\ell }(t)$ (blue) and $C/\kappa (t)$ (dashed green) are shown. Because the zonal flows suppress the turbulent particle flux, the radial density profile stops to evolve (actually, it still evolves a bit since the particle flux is not exactly zero, but extremely slowly). Indeed, $C/\kappa$ becomes almost constant at the end of the simulation. Hence, we end up with a state that is ‘frozen’ by the zonal flows. Note that, somewhat counter-intuitively, because of the dependence of the adiabaticity parameter to temperature through parallel resistivity, a higher physical temperature makes the system more easily frozen by zonal flows.

To quantify the level of zonal flows, we can measure the fraction of zonal kinetic energy $\varXi _{\mathcal{K}}\equiv \mathcal{K}_{ZF}/\mathcal{K}$ , where $\mathcal{K}_{ZF}$ is the kinetic energy of zonal modes and $\mathcal{K}$ is the total kinetic energy of the system. The fraction, shown in figure 12(b), starts from low values at the beginning of the simulation and reaches almost $70\,\%$ at the end, even though it still continues to increase.

In some simulations with a sufficiently large 2-D box, we observe that zonal flow formation can be localised near the edges of the system, while the centre is still turbulent. Then, over time, zonal flows progressively cover the domain and the turbulence becomes completely suppressed. We leave the characterisation of this ‘freezing front’ and its comparison with that of turbulence, as well as the attempt to reproduce it using the 1-D reduced model featuring an equation on the zonal velocity/shear, to future studies. However, note that the phenomenon is dependent on the value of $C$ and requires tailoring the initial density profile to be observed.