Understanding the turbulent mechanisms setting the density decay length in the tokamak scrape-off layer

Abstract

Mechanisms setting the density decay in the scrape-off layer (SOL) at the outer midplane of a tokamak plasma are disentangled using two-fluid numerical simulations in a double-null magnetic configuration and analytical estimates. Typical experimental observations are retrieved, in particular increasing intermittency of the turbulence going from the near to the far SOL, which is reflected in two different density decay lengths. The decay length of the near SOL is well described as the result of transport driven by a nonlinearly saturated ballooning instability, while in the far SOL, the density decay length is described using a model of intermittent transport mediated by blobs. The analytical estimates of the decay lengths agree well with the simulation results and typical experimental values and can therefore be used to guide tokamak design and operation.

By determining the plasma fuelling, power exhaust, impurities and neutral dynamics, the scrape-off layer (SOL) is critical to the performance of fusion devices. The SOL is the outermost plasma region of a tokamak, characterised by magnetic field lines that intersect the wall, bounded on the inner side by the last closed flux surface (LCFS) and on the outer side by the vessel wall. In the SOL, the turbulent fluctuations are of order unity, meaning there is no separation between equilibrium and fluctuation quantities and no significant separation in their length scales, presenting significant challenges to both simulation and analytical progress. The properties of the turbulence change radially across the SOL in both tokamaks and stellarators, as evidenced for example by an increase in the relative size of fluctuations (Kube et al. 2018; Niemann et al. 2020). In the near SOL, fluctuations are not intermittent, whereas in the far SOL their distribution has greater skewness and kurtosis, indicating intermittency (Boedo et al. 2003, 2014; Kuang et al. 2019). This is due to the existence of blobs – high density structures elongated along the magnetic field lines that propagate outwards due to their associated electric field (D’Ippolito, Myra & Zweben 2011). The properties of blobs have been measured extensively in tokamaks (Walkden et al. 2016; Zweben et al. 2016; Tsui et al. 2018), stellarators (Sánchez et al. 2003), reversed field pinches (Spolaore et al. 2004) and basic plasma experiments (Antar et al. 2001; Carter 2006; Furno et al. 2008a), however, the generation of blobs and prediction of blob-mediated transport remain open questions.

Recently, significant progress has been made in gyrokinetic (Dorf et al. 2012; Chang et al. 2017; Shi et al. 2017, 2019; Pan et al. 2018) and fluid (Dudson et al. 2009; Tamain et al. 2016; Häcker et al. 2018; Paruta et al. 2018; Stegmeir et al. 2018; Zhu, Francisquez & Rogers 2018; Riva et al. 2019) modelling of the SOL. Differences in the near and far SOL have been observed (Halpern & Ricci 2017) and blobs have been detected and tracked (Nespoli et al. 2017a; Paruta et al. 2019). This represents a step forward, building upon earlier single blob studies (Myra, Russell & D’Ippolito 2006) and multiblob models (Russell, Myra & D’Ippolito 2007; Militello & Omotani 2016; Walkden et al. 2017). Leveraging these achievements, in the present Letter we disentangle the different turbulent mechanisms in the SOL, in particular, the nature of the fluctuations in the near and far SOL, the properties of blobs including their typical size, velocity and generation rate and the parallel transport. Ultimately, this allows us to develop a predictive model for the SOL turbulent transport and density decay lengths. These are key elements towards predicting the heat flux scale length, which is among the most critical issues for the operation of ITER and the design of DEMO (Loarte et al. 2007; Zohm et al. 2013; Donné & Morris 2018), and determining the wall recycling, impurity influx and wall erosion. We focus on the tokamak double-null (DN) magnetic configuration. Besides being of interest for DEMO (Wenninger et al. 2016), this configuration facilitates the development of simple analytical estimates because the high and low field sides (HFS/LFS) are topologically separated.

Our analytical investigation is based on the results of two-fluid, three-dimensional (3-D) numerical simulations of the SOL dynamics carried out using the GBS code (Ricci et al. 2012; Halpern et al. 2016; Paruta et al. 2018). GBS solves the drift reduced Braginskii equations (Zeiler, Drake & Rogers 1997) to evolve self-consistently the full profiles of the density, electrical potential, electron temperature and ion and electron parallel velocities with no separation between equilibrium and fluctuations. The simulation domain covers the full toroidal and poloidal angle regime and extends radially from ${\sim}17\unicode[STIX]{x1D70C}_{s}$ inside the LCFS up to the wall. A source of heat and density within the LCFS mimics heat and plasma outflow from the core. The plasma flows along the field lines while being radially transported due to turbulence until it reaches the wall, which acts as a sink.

In the cold ion, electrostatic limit, the model equations are

(0.1)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}n}{\text{d}t}=\frac{2}{B}[C(p_{e})-nC(\unicode[STIX]{x1D719})]-\unicode[STIX]{x1D735}_{\Vert }(nv_{\Vert e})+S_{n} & \displaystyle\end{eqnarray}$$

(0.2)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}t}=\frac{2B}{n}C(p_{e})-v_{\Vert i}\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D714}+\frac{B^{2}}{n}\unicode[STIX]{x1D735}_{\Vert }j_{\Vert } & \displaystyle\end{eqnarray}$$

(0.3)$$\begin{eqnarray}\displaystyle \frac{\text{d}T_{e}}{\text{d}t} & = & \displaystyle \frac{4}{3B}\left[\frac{7}{2}T_{e}C(T_{e})+\frac{T_{e}^{2}}{n}C(n)-T_{e}C(\unicode[STIX]{x1D719})\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2}{3}\frac{T_{e}}{n}0.71\unicode[STIX]{x1D735}_{\Vert }j_{\Vert }-\frac{2}{3}T_{e}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert e}-v_{\Vert e}\unicode[STIX]{x1D735}_{\Vert }T_{e}+S_{Te}\end{eqnarray}$$

(0.4)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}v_{\Vert e}}{\text{d}t}=-v_{\Vert e}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert e}+\unicode[STIX]{x1D707}\left(\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}-\frac{T_{e}}{n}\unicode[STIX]{x1D735}_{\Vert }n-1.71\unicode[STIX]{x1D735}_{\Vert }T_{e}+\unicode[STIX]{x1D708}j_{\Vert }\right) & \displaystyle\end{eqnarray}$$

(0.5)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}v_{\Vert i}}{\text{d}t}=-v_{\Vert i}\unicode[STIX]{x1D735}_{\Vert }v_{\Vert i}-\frac{1}{n}\unicode[STIX]{x1D735}_{\Vert }p_{e}, & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ the vorticity and $d_{t}\,f=\unicode[STIX]{x2202}_{t}f+\unicode[STIX]{x1D70C}_{\ast }^{-1}[\unicode[STIX]{x1D719},f]$ the convective derivative of field $f$. The differential operators are given by $\unicode[STIX]{x1D735}_{\Vert }f=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f$, $\unicode[STIX]{x1D6FB}_{\Vert }^{2}f=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f)$, $[\unicode[STIX]{x1D719},f]=\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\times \unicode[STIX]{x1D735}f)$, $C(f)=B[\unicode[STIX]{x1D735}\times (\boldsymbol{b}/B)]\boldsymbol{\cdot }\unicode[STIX]{x1D735}f/2$ and $\unicode[STIX]{x1D6FB}_{\bot }^{2}f=\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\boldsymbol{b}\times \unicode[STIX]{x1D735}f)\times \boldsymbol{b}]$ where $\boldsymbol{b}=\boldsymbol{B}/B$ is the unit vector of the magnetic field and $B$ is its norm. Unless specified otherwise, we normalise $n$, $T_{e}$, $\unicode[STIX]{x1D719}$ and $v_{\Vert e,i}$, $B$ and $t$ to $n_{0}$, $T_{e0}$, $T_{e0}/e$, $c_{s0}$, $B_{0}$ and $R_{0}/c_{s0}$, respectively, where $n_{0}$, $T_{e0}$ and $c_{s0}=\sqrt{T_{e0}/m_{i}}$ are the reference density, electron temperature and sound speed, $B_{0}$ and $R_{0}$ are the magnetic field strength and major radius at the magnetic axis. Perpendicular lengths are normalised to the ion sonic Larmor radius $\unicode[STIX]{x1D70C}_{s0}=c_{s0}/\unicode[STIX]{x1D6FA}_{ci}$, with $\unicode[STIX]{x1D6FA}_{ci}=eB_{0}/(cm_{i})$, and parallel lengths normalised to $R_{0}$. The normalised resistivity is defined based on the Spitzer resistivity as $\unicode[STIX]{x1D708}=m_{e}R_{0}/(1.96m_{i}c_{s0}\unicode[STIX]{x1D70F}_{e})$ where $\unicode[STIX]{x1D70F}_{e}$ is the electron collision time, assumed constant across the domain. We define $\unicode[STIX]{x1D707}=m_{i}/m_{e}$ and $\unicode[STIX]{x1D70C}_{\ast }=\unicode[STIX]{x1D70C}_{s0}/R_{0}$. The electron pressure is denoted $p_{e}=nT_{e}$ and dimensionless current $j_{\Vert }=n(v_{\Vert i}-v_{\Vert e})$. The coordinates $(x,y,z)$ refer to the radial and poloidal directions in units of $\unicode[STIX]{x1D70C}_{s0}$ and the toroidal angle in radians respectively. The source terms, $S_{n}$ and $S_{T_{e}}$, are Gaussian centred at a distance $11\unicode[STIX]{x1D70C}_{s0}$ inside the LCFS with half-width half-maximum (HWHM) $1.5\unicode[STIX]{x1D70C}_{s0}$ and amplitude $1.35n_{0}$. Magnetic presheath boundary conditions (Loizu et al. 2012) are used at the wall: $v_{\Vert i}=\pm \sqrt{T_{e}},v_{\Vert e}=\pm \sqrt{T_{e}}\exp (\unicode[STIX]{x1D706}-\unicode[STIX]{x1D719}/T_{e}),\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D719}=\mp \sqrt{T_{e}}\unicode[STIX]{x2202}_{x}v_{\Vert i},\unicode[STIX]{x2202}_{x}n=\mp n/\sqrt{T_{e}}\unicode[STIX]{x2202}_{x}v_{\Vert i},\unicode[STIX]{x1D714}=-(\unicode[STIX]{x2202}_{x}v_{\Vert i})^{2}\mp \sqrt{T_{e}},\unicode[STIX]{x2202}_{xx}^{2}v_{\Vert i},\unicode[STIX]{x2202}_{x}T_{e}=0$. The plus and minus signs refer to field lines entering and leaving the wall and $\unicode[STIX]{x1D706}\approx 3$. At the inner radial boundary, we use an ad hoc set of boundary conditions: $\unicode[STIX]{x2202}_{x}\,f=0$ for all fields $f$, except for $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D719}$, for which we impose $\unicode[STIX]{x1D714}=0$ and $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D706}T_{e}$. The ad hoc inner boundary conditions have no effect in the region outside the source where our analysis is performed.

The axisymmetric magnetic field is based on three infinitely long wires aligned vertically, with the current in the central wire mimicking the plasma current. The upper and lower wires carry a current ten times stronger and are located at a distance $2a$ from the central wire with $a$ the radius at the wall. We then apply a radial transformation $x\rightarrow x-x_{0}$, where $x_{0}=0.9a$, to the flux function to obtain a configuration sufficiently circular to fit in the domain. The separatrix is shown in figure 1. Our analysis is independent of the details of the magnetic geometry.

By tuning the poloidal field strength, simulations are run with local safety factor $q=(a/\unicode[STIX]{x1D70C}_{\ast })\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}z/\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}y=4.3$, 6.5, 8.6 at the LCFS outer midplane. We also scan the parallel resistivity, $\unicode[STIX]{x1D708}=1,0.1$ and 0.01, typical experimental values for the normalised resistivity. Scans over these parameters are chosen since they are the most important controls on the plasma dynamics. We fix $\unicode[STIX]{x1D70C}_{\ast }^{-1}=500$ and $\unicode[STIX]{x1D707}=200$. The domain size is $L_{x}=120$ and $L_{y}=2\unicode[STIX]{x03C0}a=800$. Starting from a uniform initial state, the 3-D fields are evolved in time until a steady state is reached where the plasma influx from the source is balanced by losses to the wall and the toroidally averaged fields fluctuate around a constant value. We perform our analysis on this quasi-steady state. Details of the numerical implementation can be found in Paruta et al. (2018).

A snapshot of $n$ and $\bar{n}$ on a poloidal plane for the $\unicode[STIX]{x1D708}=1$, $q=6.5$ simulation is shown in figure 1 (we use an overbar and a tilde to denote the time and toroidal average and the fluctuating component of all quantities, e.g. $n=\bar{n}+{\tilde{n}}$). We observe the turbulent LFS and quiescent HFS, experimentally observed in DN configurations (LaBombard et al. 2016). Removing the interchange instability drive (the curvature term in (0.2)) drastically reduces turbulence on the LFS, indicating that ballooning instabilities are the primary driver of turbulence in this region. Removing the Kelvin–Helmholtz drive (the $[\unicode[STIX]{x1D719},\unicode[STIX]{x1D714}]$ term in the convective derivative in (0.2)) suppresses the turbulence on the HFS but has a small impact on the LFS, indicating velocity shear and the Kelvin–Helmholtz instability play a dominant role only in the HFS.

Since most of the heat is exhausted on the LFS, we focus this Letter on the mechanisms determining the density decay length at the outer midplane. The density decay in our simulations cannot be properly described by a single exponential decrease, rather, it can be fitted with two exponentials characterised by a shorter decay length $L_{n}$ near the LCFS and longer decay length $L_{n}^{\prime }$ in the far SOL (see figure 2). Such a double decay length is a typical observation in a DN configuration, e.g. on C-Mod (LaBombard et al. 2016) and MAST (Riva et al. 2019), as well as in single null (SN) (Carralero et al. 2017; Kuang et al. 2019) and limited configurations (Horacek et al. 2016). It has also been observed in two-fluid simulations (Francisquez, Zhu & Rogers 2017). The difference in the scale length is reflected in different turbulent properties in the near and the far SOL. As observed experimentally (Boedo et al. 2003; D’Ippolito et al. 2011; Kuang et al. 2019), the fluctuation distribution is close to Gaussian in the near SOL with increasing skewness and kurtosis, indicative of intermittency, in the far SOL (figure 2). In the following, we identify the two different mechanisms setting $L_{n}$ and $L_{n}^{\prime }$. We call the width of the inner SOL $\unicode[STIX]{x1D6E5}$, this is the distance over which the density decays steeply. We refer to the density, temperature and radial turbulent particle flux at the separatrix by $\bar{n}$, $\bar{T_{e}}$ and $\unicode[STIX]{x1D6E4}$, and at the entrance of the far SOL (a distance $\unicode[STIX]{x1D6E5}$ from the separatrix) by $\bar{n}^{\prime }$, $\bar{T_{e}}^{\prime }$ and $\unicode[STIX]{x1D6E4}^{\prime }$, both at the outer midplane.

Figure 1. A snapshot of $n$ (a) and $\bar{n}$ (b) on a poloidal plane for the $\unicode[STIX]{x1D708}=1$, $q=6.5$ simulation. The separatrix is shown in black. A comparison of the simulation results and model prediction is shown in (c).

Figure 2. Radial profile of $\bar{n}$ at the outer midplane with the fitted exponentials (a), profiles of the skewness and kurtosis (b) and total turbulent transport and transport due to blobs (c). Above are two PDFs of the normalised density fluctuation, ${\tilde{n}}/\unicode[STIX]{x1D70E}_{n}$, evaluated in the corresponding radius range, where $\unicode[STIX]{x1D70E}_{n}$ is the standard deviation evaluated locally.

We start by looking at the near SOL. Since the turbulent radial flux in the near SOL is not intermittent, we estimate the flux based on the development and saturation of a linear instability driven by a background radial gradient in density and temperature. We then match the predicted flux, $\unicode[STIX]{x1D6E4}$, to the turbulent flux across the LCFS, $\unicode[STIX]{x1D6E4}_{\text{LCFS}}$, to find $L_{n}$. It should be noted that although ${\tilde{n}}(k_{y})$ has a broad spectrum, $\unicode[STIX]{x1D6E4}(k_{y})$ has a clear peak (Podestà et al. 2008).

The turbulent flux can be written $\unicode[STIX]{x1D6E4}=\langle \overline{{\tilde{n}}\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}}/B\rangle _{y}$, where the poloidal, $y$, average is evaluated over $45^{\circ }$ centred around the outer midplane. The density fluctuation can be estimated by noticing that linear instabilities saturate when the gradient of the fluctuations becomes comparable to the background density gradient, hence locally removing the turbulence drive (Ricci, Rogers & Brunner 2008; Ricci & Rogers 2013), that is $\unicode[STIX]{x2202}_{x}{\tilde{n}}\sim \unicode[STIX]{x2202}_{x}\bar{n}$ or equivalently $k_{x}{\tilde{n}}\sim \bar{n}/L_{n}$, with $k_{x}$ the typical radial wavenumber of the perturbation, in agreement with ${\tilde{n}}/\bar{n}$ in the simulations. We relate $\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}$ to ${\tilde{n}}$ by balancing the leading-order terms in the continuity equation, equation (0.1): $\unicode[STIX]{x1D6FE}{\tilde{n}}=\unicode[STIX]{x1D70C}_{\ast }^{-1}\unicode[STIX]{x2202}_{y}\tilde{\unicode[STIX]{x1D719}}\unicode[STIX]{x2202}_{x}\bar{n}/B$ where $\unicode[STIX]{x1D6FE}$ is the linear growth rate of the instability driving the transport. The simulation test mentioned above shows that ballooning modes drive turbulence on the LFS. For these modes, as well as for drift waves (Rogers & Dorland 2005), non-local linear theory shows that $k_{x}=\sqrt{k_{y}/L_{n}}$ (Ricci & Rogers 2013). By considering the linear instability that maximises the transport, the turbulent flux $\bar{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D70C}_{\ast }\bar{n}(\unicode[STIX]{x1D6FE}/k_{y})_{\text{max}}$ follows.

In order to evaluate $(\unicode[STIX]{x1D6FE}/k_{y})_{\max }$, we linearise (0.1)–(0.5) assuming $C(f)\sim \unicode[STIX]{x2202}_{y}$, since the flux surfaces are approximately vertical in most of the region we are considering, and neglecting radial variation of the perturbation and poloidal variation of the equilibrium. We take $k_{\Vert }=2/q$, where 2 is the minimum parallel mode number, expected from ballooning stability and observed in the simulations. The linearised system that we obtain corresponds to that of the simple magnetic torus geometry (Poli et al. 2008). We take $\unicode[STIX]{x1D702}=L_{n}/L_{Te}=0.77$, the theoretically expected value (Ricci et al. 2008), which is similar to the simulations. Using $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{\text{LCFS}}$ we find numerically the $k_{y}$ and $L_{n}$ shown in figure 3 for three values of $\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. The estimates of $L_{n}$ correspond well to the simulations, as shown in figure 1.

To understand the decrease in $k_{y}$ and increase in $L_{n}$ with increasing $\unicode[STIX]{x1D708}$ and $q$ (which has been observed experimentally Nespoli et al. 2017b), we consider the limit (valid for typical parameters) in which the resistive ballooning mode is dominant: $R/L_{n}\gg 1$, $\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FE}\ll \unicode[STIX]{x1D708}$ and, to avoid coupling with sound waves and drift waves, $k_{\Vert }\ll \unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D714}_{\ast }\ll \unicode[STIX]{x1D6FE}$ where $\unicode[STIX]{x1D714}_{\ast }=k_{y}R/L_{n}$. Our dispersion relation reduces to $\unicode[STIX]{x1D6FE}^{2}-\unicode[STIX]{x1D6FE}_{i}^{2}+\unicode[STIX]{x1D6FE}k_{\Vert }^{2}/(\unicode[STIX]{x1D708}k_{y}^{2})=0$, where $\unicode[STIX]{x1D6FE}_{i}=\sqrt{2R(1+\unicode[STIX]{x1D702})/L_{n}}$ is the ideal ballooning growth rate, capturing well the strong transport limit (Halpern et al. 2014). Expressed in physical units to make explicit the $R_{0}$ dependence,

(0.6)$$\begin{eqnarray}L_{n}=\frac{2}{3}(1+\unicode[STIX]{x1D702})\bar{c_{s}}^{2}\left(\frac{q\bar{n}}{2\unicode[STIX]{x1D6E4}_{\text{LCFS}}}\right)^{4/3}\left(\frac{m_{e}}{m_{i}1.96\unicode[STIX]{x1D70F}_{e}}\right)^{2/3}\unicode[STIX]{x1D6FA}_{ci}^{-4/3}R_{0}^{1/3},\end{eqnarray}$$

which is of the order of 1 mm for typical experimental parameters in TCV and C-Mod (LaBombard et al. 2016; Vianello et al. 2019). We can use this equation to estimate $k_{y}$ by noticing that the peak $\unicode[STIX]{x1D6FE}/k_{y}$ occurs approximately where the damping term and ballooning drive are equal. Using $\unicode[STIX]{x1D6FE}\sim \unicode[STIX]{x1D6FE}_{i}$, we find $k_{y}=(2^{2/3}/3^{1/4})[\bar{n}\unicode[STIX]{x1D6FA}_{ci}^{2}(1.96m_{i}\unicode[STIX]{x1D70F}_{e})/(q^{2}\unicode[STIX]{x1D6E4}_{\text{LCFS}}R_{0}^{2}m_{e})]^{1/3}$ in physical units, the full numerical result is shown in figure 3.

Figure 3. Predicted $L_{n}$, $k_{y}$, $v_{b}$, $\unicode[STIX]{x1D6E4}^{\prime }$ and $L_{n}^{\prime }$ as a function of $\unicode[STIX]{x1D708}$, $q$ and $\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. All reference quantities are taken at the LCFS.

We now turn to the far SOL, where the fluctuation distribution is heavy tailed, indicating intermittent turbulence and indeed observation of the simulation results reveals the presence of coherent structures of high plasma density (blobs) that propagate outwards due to their self-generated $\boldsymbol{E}\times \boldsymbol{B}$ velocity.

We use a pattern recognition algorithm described in Paruta et al. (2019) to track the blobs (defined here as coherently propagating structures of amplitude greater than 2.5 times the standard deviation of $n$) and measure their size, amplitude and velocity. Following Nespoli et al. (2017a), we calculate the fraction of the cross-field transport due to blobs by assuming a 2-D Gaussian density distribution of each blob in the poloidal plane with a peak density fluctuation $n_{b,i}$ and radial and poloidal HWHM $a_{x,i}$ and $a_{y,i}$, where $i$ is the blob index. The blob flux is calculated by $\unicode[STIX]{x1D6E4}_{b}(x,y)=\sum _{i}n_{b,i}v_{b,i}\exp [(x-x_{b,i})^{2}/(2a_{x,i}^{2})+(y-y_{b,i})^{2}/(2a_{y,i}^{2})]$, where the sum is carried out over all blobs and $(x_{b,i},y_{b,i})$ are the blobs’ centre of mass. We find that blob transport dominates in the far SOL (figure 2), consistent with the result of Nespoli et al. (2017a) and previous experimental works that found blobs to contribute an order unity fraction of the particle flux (Boedo et al. 2003, 2014; D’Ippolito et al. 2011).

We now predict the flux due to blobs using only the near SOL properties. For this purpose, we express the blob flux averaged in the poloidal plane (Russell et al. 2007) $\unicode[STIX]{x1D6E4}^{\prime }=\langle \unicode[STIX]{x1D6E4}_{b}\rangle _{x,y}=\unicode[STIX]{x1D70E}_{b}f_{b}v_{b}$ in terms of $\unicode[STIX]{x1D70E}_{b}$ the average density inside a blob, $f_{b}$ the blob packing fraction (ratio of area covered by blobs to total SOL area) and $v_{b}$ the average blob velocity. The $x$ average is taken from the LCFS to the wall. We address each of these quantities in turn.

We estimate $\unicode[STIX]{x1D70E}_{b}=2n_{b}/\text{ln}(2)$, i.e. as the ratio of the average number of particles in a blob, $2\unicode[STIX]{x03C0}n_{b}a_{x}a_{y}/\text{ln}(2)$, where we assume blobs have on average a Gaussian shape with HWHM $a_{x}$ and $a_{y}$ and peak density $\bar{n}^{\prime }+n_{b}$, and the average blob area, $A_{b}=\unicode[STIX]{x03C0}a_{x}a_{y}$. Since $\bar{n}^{\prime }$ decreases radially, $n_{b}/\bar{n}^{\prime }$ remains approximately constant over the blob lifetime despite parallel draining, so we combine the definition of a blob and the estimate of the ${\tilde{n}}/\bar{n}$ in the near SOL to estimate $n_{b}\sim 3{\tilde{n}}\sim 3\bar{n}/(L_{n}k_{x})$.

The packing fraction $f_{b}=N_{b}A_{b}/(A_{SOL})$, where $N_{b}$ is the number of blobs, requires an estimate of the blob size. We observe that the blob size remains approximately constant as the blobs propagate (as observed experimentally in AUG Carralero et al. 2015) and that blobs tend to be circular, maximising their Kelvin–Helmholtz stability (Ricci & Rogers 2013), so we estimate their size as the geometric mean of the near SOL eddy dimensions $a_{x}\sim a_{y}\sim \unicode[STIX]{x03C0}/(2\sqrt{k_{x}k_{y}})$. We infer, therefore, from the results for the near SOL analysis that blob size increases with resistivity, a trend observed both in our simulations (figure 4) and experimentally Vianello et al. (2019).

Figure 4. Distribution of the vertical size (a) and radial velocity (b) of the blobs at $\unicode[STIX]{x1D708}=1$ and $\unicode[STIX]{x1D708}=0.01$ with $q=6.5$. The mean sizes and velocities are shown with a solid line and the predictions with a dashed line.

We now turn to the estimation of $N_{b}$. In steady state, the blob generation and loss rates are equal. Since blobs are generated from instabilities of wavelength $2\unicode[STIX]{x03C0}/k_{y}$, we expect the generation rate to be proportional to $L_{y}k_{y}/(2\unicode[STIX]{x03C0})$. The generation time scale has previously been proposed as determined by poloidal flow shear (Furno et al. 2008b) or a combination of flow shear and mode phase velocity (D’Ippolito et al. 2011; Fuchert et al. 2016). In the presence of hot ions, strong $\boldsymbol{E}\times \boldsymbol{B}$ flow shear may be present (Zhu, Francisquez & Rogers 2017; Paruta et al. 2018). However, in our simulations, we find a flow shear time scale, $\unicode[STIX]{x2202}_{x}v_{E\times B}$, almost an order of magnitude less than the observed generation time scale and not to scale with $\unicode[STIX]{x1D708}$ and $q$, as observed. We reason that blobs are created because the linear instability saturates as the local density gradient is removed and the resulting density perturbation moves outwards without the streamer being refilled from the core, a case which was studied in a basic plasma physics device in Müller et al. (2009). Hence, the generation rate is limited by the time taken for the blob to travel one radial wavelength of the driving instability, $4a_{x}/v_{b}$, allowing for the density gradient to be re-established, which is consistent with the simulations. Taking the blob lifetime as the time taken to cross the domain, the loss rate is $N_{b}v_{b}/L_{x}$. Hence, $N_{b}=4\unicode[STIX]{x03C0}^{2}L_{x}L_{y}/(k_{x}k_{y})$. Using the above relations, we find $f_{b}\approx \unicode[STIX]{x03C0}/16$, independent of the SOL parameters. While the universality of $f_{b}$ is well supported by the simulation data, the predicted value is an overestimate (likely because we assume all blobs cross the entire radial domain) and a better estimate is $f_{b}=0.1$.

The average blob velocity, $v_{b}$, is deduced from the average blob size, $a_{x}$ and $a_{y}$, according to the well-studied size–velocity scaling relations derived using the two-region model (Myra et al. 2006; D’Ippolito et al. 2011; Tsui et al. 2018; Paruta et al. 2019) (figure 4). The normalised blob velocity, $\hat{v}=\text{Im}(\hat{\unicode[STIX]{x1D714}})\hat{a}^{1/2}$, depends on the normalised frequency $\hat{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FE}_{b}$ (where $\unicode[STIX]{x1D6FE}_{b}^{2}=2\bar{T_{e}}^{\prime }\unicode[STIX]{x1D70C}_{\ast }^{-1}n_{b}/[a_{x}\bar{n}^{\prime }]$ represents the local ballooning drive), which is determined by the dispersion relation $1+\hat{\unicode[STIX]{x1D714}}^{2}+\text{i}\hat{\unicode[STIX]{x1D714}}\unicode[STIX]{x1D6E9}/\unicode[STIX]{x1D6EC}=0$, for $\unicode[STIX]{x1D6EC}>1$, or $1+(1+f^{2})\hat{\unicode[STIX]{x1D714}}^{2}+\text{i}\unicode[STIX]{x1D6E9}\hat{\unicode[STIX]{x1D714}}=0$, for $\unicode[STIX]{x1D6EC}<1$, where $\unicode[STIX]{x1D6EC}=\unicode[STIX]{x1D708}\bar{n}^{\prime }L_{\Vert 1}^{2}/[L_{\Vert 2}\unicode[STIX]{x1D70C}_{s}]>1$ is the resistivity parameter, $L_{\Vert 1}$ is the field line length from the midplane to the region of maximum flux fanning, $L_{\Vert 2}$ is the field line length from the region of maximum flux fanning to the wall (Paruta et al. 2019) and $f$ is the flux tube fanning (Myra et al. 2006). The blob size parameter, $\unicode[STIX]{x1D6E9}=\hat{a}^{5/2}$, where $\hat{a}=a_{b}/a_{\ast }$ is the normalised blob size, with $a_{b}=(2a_{y}/\unicode[STIX]{x03C0})^{4/5}a_{x}^{1/5}$ and $a_{\ast }=[2\unicode[STIX]{x1D70C}_{s}^{4}L_{\Vert 2}^{2}n_{b}/(a_{x}\bar{n}^{\prime }\unicode[STIX]{x1D70C}_{\ast })]^{1/5}$. Finally the blob velocity is given by $v_{b}=0.5\hat{v}v_{\ast }$, with the reference velocity $v_{\ast }=\unicode[STIX]{x1D70C}_{s}[2\unicode[STIX]{x03C0}^{2}a_{x}^{2}\unicode[STIX]{x1D70C}_{s}^{2}\unicode[STIX]{x1D70C}_{\ast }^{2}L_{\Vert 2}n_{b}/(\bar{n}^{\prime }a_{y}^{2})]^{1/5}$ and the factor 0.5, obtained by comparing the scaling with the simulation results, accounting for the fact that our estimate is an upper limit that neglects various mechanisms slowing the blobs (Tsui et al. 2018).

Finally, we determine $L_{n}^{\prime }$ by balancing the divergence of the blob flux with the divergence of the parallel flow $L_{n}^{\prime }=\unicode[STIX]{x1D6E4}^{\prime }L_{\Vert }^{\prime }/(\bar{n}^{\prime }\bar{c_{s}}^{\prime })$ with $L_{\Vert }=L_{\Vert 1}+L_{\Vert 2}$, $\bar{n}^{\prime }=\bar{n}\exp (-\unicode[STIX]{x1D6E5}/L_{n})$ and $\bar{c_{s}}^{\prime }=\bar{c_{s}}\exp (-\unicode[STIX]{x0394}\unicode[STIX]{x1D702}/[2L_{n}])$. We remark that $L_{n}^{\prime }$ depends only weakly on $\unicode[STIX]{x1D6E5}$ since $\unicode[STIX]{x1D6E4}^{\prime }$ also scales approximately with $\bar{n}^{\prime }\bar{c_{s}}^{\prime }$. For typical experimental parameters, most blobs are in the $\unicode[STIX]{x1D6EC}>\unicode[STIX]{x1D6E9}$ regime, for which

(0.7)$$\begin{eqnarray}L_{n}^{\prime }=\frac{7.3f_{b}L_{\Vert }\unicode[STIX]{x1D6FA}_{ci}^{0.0167}(R_{0}m_{i}\unicode[STIX]{x1D70F}_{e})^{0.00833}\unicode[STIX]{x1D6E4}_{\text{LCFS}}^{1.04}}{m_{e}^{0.00833}\bar{c_{s}}^{1.05}\bar{n}^{1.04}(1+\unicode[STIX]{x1D702})q^{0.0167}\unicode[STIX]{x1D70C}_{s}^{\prime 0.1}}\end{eqnarray}$$

is of the order of several mm for typical experimental parameters. We note that, since $\unicode[STIX]{x1D6EC}\propto R_{0}$ and $\unicode[STIX]{x1D6E9}\propto R_{0}^{1.46}$, larger device will likely have blobs in the $\unicode[STIX]{x1D6E9}>\unicode[STIX]{x1D6EC}$, $\unicode[STIX]{x1D6EC}>1$ regime, for which

(0.8)$$\begin{eqnarray}L_{n}^{\prime }=\frac{94.8f_{b}L_{\Vert 1}^{2}L_{\Vert }\unicode[STIX]{x1D6E4}_{\text{LCFS}}^{5/3}m_{i}^{17/15}\unicode[STIX]{x1D70F}_{e}^{2/15}\unicode[STIX]{x1D70C}_{s}^{\prime 7/5}\unicode[STIX]{x1D6FA}_{ci}^{19/15}}{R_{0}^{28/15}\bar{n}^{5/3}m_{e}^{17/15}c_{s}^{14/5}\unicode[STIX]{x1D70C}_{s}^{2/5}(1+\unicode[STIX]{x1D702})^{7/5}q^{34/15}}.\end{eqnarray}$$

In figure 3 we show $v_{b}$, $\unicode[STIX]{x1D6E4}^{\prime }$ and $L_{n}^{\prime }$ as a function of $\unicode[STIX]{x1D708}$, $q$ and $\unicode[STIX]{x1D6E4}_{\text{LCFS}}$. We observe that $\unicode[STIX]{x1D6E4}^{\prime }$ increases with $\unicode[STIX]{x1D708}$ and with $q$, primarily due to variation in $v_{b}$ and to a lesser extent $\unicode[STIX]{x1D70E}_{b}$, as suggested in Russell et al. (2007). The increase in $v_{b}$ follows from the $\hat{a}-\hat{v}$ scaling. Such an increase has also been observed in gyrofluid simulations (Häcker et al. 2018). The increase in $L_{n}^{\prime }$ with resistivity is well documented experimentally (D’Ippolito et al. 2011). The predicted $L_{n}^{\prime }$ is compared to the simulation result in figure 1. As for the near SOL, we find good agreement between theory and simulation.