Dynamical simulations of boundary plasma turbulence in divertor geometry

Direct comparisons between numerical simulations and the measured plasma fluctuations and transport are presented by performing nonlinear two-fluid simulations with the BOUT code (Xu X Q and Cohen R H 1998 Contrib. Plasma Phys. 38 158). BOUT models boundary-plasma turbulence in a realistic divertor geometry using modified Braginskii equations for plasma vorticity, density (ni), electron and ion temperature (Te, Ti) and parallel momenta. The BOUT code solves for the plasma fluid equations in a 3D toroidal segment, including the region somewhat inside the separatrix and extending into the scrape-off layer; the private flux region is also included. In this paper, the physics of resistive X-point turbulence and its relation to flow shear generation is discussed. We present comparisons between the boundary plasma turbulence observed in the BOUT code and experiments on DIII-D (Luxon J L et al 1986 Int. Conf. on Plasma Physics and Controlled Nuclear Fusion (Vienna: IAEA) p 159), the National Spherical Torus Experiment (Peng Y-K M 2000 Phys. Plasmas 7 1681), and C-Mod (Hutchinson I H et al 1994 Phys. Plasmas 1 1511). In an L-mode discharge in the DIII-D tokamak, both BOUT simulations and beam emission spectroscopy show a similar flow pattern and blob size across the last closed flux surface. In an L-mode discharge, both BOUT simulations and gas puff imaging show similar filament structures along the field line and similar frequency spectrum at the outboard midplane. In simulations of the quasi-coherent mode in the EDA regime of C-Mod, the particle flux measured from BOUT simulation is consistent with Langmuir probe measurements on C-Mod at the midplane near the separatrix. The qualitative comparisons thus indicate that BOUT contains much of the relevant physics for boundary plasma turbulence in the experimentally relevant X-point divertor geometry of present-day tokamaks and spherical tori.


Introduction
The performance of tokamaks and other toroidal magnetic fusion devices depends crucially on the dynamics of the boundary region, i.e. the transition region from the hot core plasma through the separatrix to the material surface of the first wall, as shown in figure 1. Plasma turbulence, and the resulting anomalous cross-field plasma transport, are physical processes in the boundary region, affecting core plasma confinement (e.g. H-mode), the density limit and plasma-wall interactions [1]. The plasma boundary region has a number of physics attributes which make it distinct from the core: relatively low temperature, large radial gradients, high neutral-gas and impurity densities, proximity of open and closed flux surfaces, presence of X-point and sheath physics in the scrape-off layer (SOL). The large radial gradients tend to drive turbulent fluctuations which are a larger percentage of background values than in the core plasma.
Strong boundary turbulence has been observed in nearly all magnetic confinement devices. The turbulence has many common features, and a great deal of experimental data have been obtained over the past 20 years on, for example, fluctuation levels, spectra, correlation lengths and scalings [2], but until recently these data could not be understood from first principles. With the recent development of three-dimensional nonlinear codes, such as BOUT, it has become possible to make a direct computation of boundary turbulence, and validating these codes with experiments has now begun [3]- [5].
Predictive simulation of boundary turbulence from fundamental physics models is an important, unique, and daunting challenge owing to the special properties of the boundary plasma, its importance to an overall understanding of fusion plasmas, and the vastness of its spatial and temporal scales. A critical first task is to demonstrate that simulations are able to reproduce the phenomena observed in real magnetic confinement devices. In this paper we present initial comparisons between the boundary plasma turbulence observed in the BOUT code and experiments on NSTX (National Spherical Torus Experiment), DIII-D and C-Mod.
Such comparisons are enabled not only by advances in simulation capability, but also by development of new experimental diagnostics that can probe the structure of the turbulence.

53.3
Two such diagnostics are beam emission spectroscopy (BES) [6] and gas puff imaging (GPI) [7]. In these diagnostics, the working principles are similar. The contrast and brightness of the light fluctuation in the plasma is increased by the introduction of either a neutral beam (BES) or a localized gas puff (GPI) in the field of view of the imaging diagnostic. In BES, the 32 spatial channels are arranged to image a 5 × 6 cm 2 (radial × poloidal) region in the plasma cross section, at a nominal 1 cm spatial resolution and separation. In GPI, a camera with a 10 µs gating time images D α (656 nm) light using an optical system with ∼2-3 mm spatial resolution over ∼6 × 6 cm 2 field. Turbulence can be seen with camera gating of ∼2 µs/frame as structure within the time-averaged envelope of the D α cloud.
The remainder of the paper is organized as follows. Section 2 discusses X-point effects on boundary turbulence and its characteristics. The simulation results for the flow shear generation are given in section 3. The initial validations with experimental measurements is discussed in section 4. Finally, a summary of this paper is presented in section 5.

Dynamical equations
BOUT models boundary-plasma turbulence in a realistic divertor geometry using modified Braginskii equations [8] for plasma vorticity , density (n i ), electron and ion temperature (T e , T i ) and parallel momenta. The dynamical equations are discussed in [3]. Here we examine the role of magnetic pumping in determining the magnitude of E r in the full (X-point) geometry; we believe that such an analysis has not been presented before. The basic equation to be discussed is equation (1) of [3], and which is rewritten as follows: The magnetic pumping term is the one with the pressure anisotropy, which is caused by flows in the plasma confined by a spatially varying magnetic field; it has been derived previously in collisional regimes [9]: Definitions of various quantities are as follows:

53.4
Here The symbol tilde represents the fluctuation quantities. Also, µ ii is the classical diffusion coefficient, and ν ei is electron collision frequency. The magnetic pumping term defined by equation (2) is important because it damps shear in the plasma flow. While the term ∇ ln n · V P i in (P − P ⊥ ) i,e is larger than other terms in equation (2) by order of k ⊥ R 0 for the perturbed quantities, its flux surface average vanishes. The magnetic pumping term makes a negligible contribution to linear instability because it is on the order of 1/k ⊥ R 0 smaller than other dominant linear curvature drives, such as the ∇P term in equation (1).
Averaging equation (1) over a flux surface yields In order to efficiently simulate turbulence with short perpendicular wavelengths compared with parallel wavelengths (i.e. for wavenumber k k ⊥ ), we choose field-line-aligned ballooning coordinates, (x, y, z), which are related to the usual flux coordinates (ψ, θ, ϕ) by the relations x = ψ − ψ s , y = θ, and z = ϕ − ν(x, y) dy. Here ν ≡ a e B t /RB θ is the local pitch, where a e is the effective minor radius, R is the major radius, and B t,θ are the toroidal and poloidal magnetic fields, respectively. The precise definition of the local pitch ν is given in section 2.2. The partial derivatives are: ∂/∂ψ = ∂/∂x − ( dy ∂ν/∂ψ)∂/∂z, ∂/∂θ = ∂/∂y −ν ∂/∂z, ∂/∂ϕ = ∂/∂z, and ∇ = (RB θ /a e B)∂/∂y. The magnetic separatrix is denoted by ψ = ψ s . Here the key ballooning assumption is |∂/∂y| |ν ∂/∂z| and ∂/∂θ −ν ∂/∂z. In this choice of coordinates, y, the poloidal angle, is also the coordinate along the field line. Using the ballooning coordinate system and keeping only the linear terms in ∇ ln n · V P i , we obtain Here the mixed coordinates for ∂/∂ψ and ∂ have been used for brevity. Substituting equation (4) into equation (3) yields The reason for keeping only the linear terms in ∇ ln n · V P i is that the nonlinear term is 1/k ⊥ R 0 smaller than the nonlinear E × B term on the left-hand side. The third term on the right-hand The flux surface average is performed over y and z coordinates. Because of the definite sign of the second term on the right-hand side, this equation clearly shows that magnetic pumping damps out flow shear. We note that a similar damping mechanism in the plateau-banana regime yields a nonlinear dependence on the poloidal velocity. This has been shown to lead to bifurcated solutions in the poloidal momentum balance equation in tokamaks when orbit loss is taken into consideration in the vicinity of separatrix [10].

53.6
The main complication in a diverted tokamak is the presence of the X-point, which significantly changes the structure of the fluctuations due to magnetic shear (field-line fanning) effects. In an X-point divertor geometry, the local ν(ψ, θ) profiles, as shown in figure 2 for a typical single-null divertor, vary along the magnetic field line as well as in the radial direction, and there exist steep parallel and radial gradients of ν(ψ, θ) near the X-point. There are two effects: (i) local magnetic shear becomes large near the X-point; (ii) the magnetic connection length becomes infinite when a flux tube passes near an X-point.
For single-null divertor geometry, the curvature is destabilizing on the low-magnetic-field side (outside torus) and stabilizing on the high-field side. The dependence of the normal curvature along a magnetic field line is almost a negative constant on the low-magnetic-field side with a negative dip near the X-point region and is almost a positive constant on the high-field side with a positive peak near the X-point region. The transition from the negative constant to a positive one occurs rather abruptly at the top of the machine. Frequently 'single-null' geometry is really unbalanced double null, and the second X-point in the far SOL region can still significantly affect the shear and ν(ψ, θ) as in figure 2.

Resistive X-point modes
Recent studies of resistive ballooning modes in the boundary plasma of diverted tokamaks have been performed within the framework of a collisional fluid model [3,11]- [13]. It is shown that the large magnetic shear and small poloidal field in the X-point region act to increase the local wavenumber, and hence the importance of the resistivity, near the X-point, as shown in figure 3. This effect can be viewed as a consequence of the fanning of the flux tubes in the X-point region [14], and has no simple analogue in circular flux geometry where the magnetic shearŝ is constant on a flux surface. As a result of the locally enhanced resistive effects, the mode can line-bend across the X-point, avoiding the good curvature region on the inside of the torus. This 'disconnection' of the eigenmodes profoundly influences the unstable spectrum. A new class of modes called resistive X-point (RX) modes exploits this synergism between resistivity and the X-point geometry, giving rise to robust growth rates at moderate-to-low mode numbers, for which resistive effects would otherwise be negligible.

Edge current and peeling-ballooning modes
In H-mode discharges, the sharp pressure gradients in the pedestal region can drive large bootstrap currents which provide an additional source of free energy to drive MHD instabilities. This current drives 'peeling' or edge-localized external kink modes [15,16] which, due to the expected sharp current gradient near the separatrix, can be unstable even at relatively high values of the toroidal mode number (n). In addition to driving peeling modes, the large bootstrap current in the edge region also lowers the local magnetic shear, and in shaped discharges this can allow second stability access to high-n ballooning modes. The stability physics is further complicated by the coupling of peeling modes to pressure-driven ballooningtype instabilities which occurs at finite n [15,17]. As a result, intermediate 3 < n < ∼ 40 coupled peeling-ballooning modes are often the limiting ideal MHD instability in the pedestal. It has been proposed that these peeling-ballooning modes are responsible for edge-localized modes (ELMs), and that they provide constraints on the height of the H-mode pedestal.  Predictions from peeling-ballooning stability calculations have yielded encouraging agreement with observed ELM onset times and penetration depth, and variation in pedestal characteristics with plasma shape [17].
Most previous studies of peeling-ballooning modes have been carried out using the ideal MHD model, employing equilibria which are cut off just inside the separatrix, and using a vacuum to model the resistive plasma in the scrape-off layer. While this likely provides a reasonable approximation for exploring the linear behaviour of modes driven predominantly inside the SOL, important additional physics such as ExB and diamagnetic stabilization, resistivity, and geometric effects due to divertor geometry may also play a role. In order to study these effects, and to begin exploring the nonlinear dynamics of ELMs, we have incorporated edge current into the BOUT physics model. In the Braginskii equations, the dominant 'kink' term due to the equilibrium parallel current J 0 enters in the vorticity equation: Bootstrap models are used to construct equilibria including parallel current J 0 , and these equilibria are used in BOUT to assess the impact of current and peeling modes. A preliminary set of simulations [18] has found that above a threshold current, fast-growing modes are observed with a linear structure that extends from the sharp-gradient region of the pedestal out to the separatrix. The modes exhibit a peeling-ballooning structure similar to that observed in linear ideal MHD analysis. In the early nonlinear phase, finger-like structures extend across the outer midplane and transport particles into the SOL. Further study of the nonlinear dynamics of these modes with BOUT is an important direction for future study.

Numerical scheme
The BOUT code solves the plasma fluid equations in a 3D toroidal segment, including the region somewhat inside the separatrix and extending into the SOL; the private flux region is also included. A finite-difference method is used, and the resulting difference equations are solved with a fully implicit Newton-Krylov solver PVODE [19,20], a parallelized ordinary differential equation (ODE) solver. BOUT is parallelized via a poloidal domain decomposition model that uses the MPI (message passage interface) system. The parallel implementation is straightforward and efficient: one or several poloidal meshes with the entire radial-toroidal plane are stored on each processor. At the end of a time step, the data in the domain boundary planes are passed to its physically neighbouring processor. Because of this parallel paradigm, the amount of message passing scales linearly with the problem size. For a typical run with 64 processors, the communication time is less than 1%. With poloidal flux, ψ, normalized to unity on the separatrix, we typically take the inner simulation boundary condition to be ψ c = 0.9-0.95 and the outer boundary at ψ w = 1.05-1.

Simulation results for flow generation
To simulate boundary plasma turbulence and validate with the corresponding experiments, the BOUT code uses realistic X-point magnetic and plasma profiles [3,12]. The background magnetic field structure is obtained from an MHD equilibrium code (usually EFIT [21]) for a typical shot. The plasma profiles are obtained by taking density and temperature as analytic (modified tanh) fits to Thomson scattering data ( figure 4). For typical DIII-D boundary plasma profiles in L-mode, the midplane values on the separatrix are T e = 60 eV, T i = 240 eV, and N i = 6.5 × 10 18 m −3 . From the given magnetic geometry and plasma profiles corresponding to a specific experimental device and discharge, the simulation is initialized with a set of small random fluctuations. The fastest growing modes dominate the initial phase of the calculation, in which the fluctuations grow at an approximately exponential rate. After this initial linear phase, the density and electrostatic potential fluctuations evolve to a saturated state, as shown in figure 5. Changes in the background density and temperature due to the turbulence are set to zero. However, the electric potential, parallel current and ion velocity are self-consistently evolved with turbulence dynamics. At t ∼ 43 µs, the unstable modes inside the separatrix enter into a nonlinear phase. After a period of adjustment, the turbulence-generated electric potential reaches a steady state. Turbulence-driven shear flows are now self-consistently included in the nonlinear BOUT boundary plasma simulations [22]. The flows are set by toroidal momentum balance (including Reynolds stress [23]) and can be considered as a nonlinear instability associated with inverse cascade of the turbulent spectra [24]- [26]. This flow was slowly damped through the magnetic pumping effect at a rate on the order of 1/k ⊥ R 0 smaller than the typical characteristic turbulence frequency, as discussed in section 2.1. A similar time history for turbulence evolution and flow generation has been observed for a shifted-circle geometry [27]. The negative E r inside the separatrix is generated by the Reynolds stress and positive E r in the far SOL is dominated by the sheath physics due to parallel particle loss. Near the separatrix in the SOL, the two mechanisms compete. This qualitative shape of E r is consistent with the speculation in an earlier L-H transition model [28]. Here the boundary condition for δφ is homogeneous Neumann, d δφ /dx = 0, at the inner core-edge boundary, and Dirichlet, δφ = 0, at the wall in the far SOL. By comparison, we find that the ion diamagnetism also contributes to the flow generation, but it is smaller than the Reynolds stress in L-mode in boundary regions where the resistive X-point mode is the dominant mode. We note the recent finding that geodesic acoustic modes (GAM) play a prominent role in flow generation in the core-edge region where the ion temperature gradient mode dominates [29]. These modes have a relatively weaker effect in the region where the present simulations are performed. In particular, GAM oscillations are not evident in figure 5(b).

53.11
outer-midplane inner-midplane λ ∼10m || λ ∼5cm ⊥ NSTX GPI View BOUT results Figure 7. Filament-like structures observed in BOUT and GPI. Images of the edge turbulence in the poloidal versus toroidal plane. Shown is the visible emission intensity for GPI obtained with 10 µs exposure and no interference filter. The high recycling layer over the centre column can be seen close to the right edge of these images (discharge 101 533). BOUT plot is the density fluctuation (reprinted by permission from [7]).

Initial validation with experimental measurements
Most properties of the turbulence from BOUT simulations in boundary plasmas are consistent with the various fluctuation measurements. From the correlation function, parameters like fluctuation levels, spectra, correlation time, poloidal correlation length and poloidal propagation velocity have been extracted from the BOUT simulations of DIII-D, NSTX and C-Mod using the diagnostic code GKV. The qualitative comparisons thus indicate that BOUT contains much of the relevant physics for the boundary plasma turbulence for the experimentally relevant X-point divertor geometry. To simplify the validation process by avoiding uncertainty associated with the models of source and sink, we concentrate on examining a more detailed description of boundary turbulence characteristics based on the fixed plasma profiles given by the averaged Thomson scattering data. Thus we will leave the important topics of the large transport events [30]- [32] and plasma profile evolution to future publications. As we discussed in figure [7], data from a discrete channel detector which viewed a 5 cm diameter area within the 60 cm diameter gas puffing port (reprinted by permission from [7]).
shear structure. Figure 6 provides a direct visualization of 'quasi-2D' structures in boundary turbulence. These structures are quasi-2D in that they are strongly elongated along the magnetic field line (see discussion below); we are looking at a slice almost perpendicular to the magnetic field. An animation is also given, showing large-scale, transient, and coherent structures (localized in time and in space) convecting upward through the observation domain at several kilometres per second and reversing outside the separatrix. The poloidal phase velocity v ph of the potential fluctuations in the midplane is in the ion drift direction and becomes negative (electron drift direction) close to the last closed flux surface (LCFS) and into the edge of the main plasma. Velocity reversal has been observed in tokamaks and stellarators [33]- [36]. The reversal of the fluctuation phase velocity with x is due to different dynamics inside the separatrix and in the SOL. The RX turbulence, mainly inside the separatrix, is driven by electron skin physics coupled with bad curvature. The RX turbulence phase velocity is in the electron drift direction. Due to the X-point physics, the RX turbulence penetrates some way into the SOL. Because the electron temperature decays faster than the ion temperature in the SOL, the ion diamagnetic drift velocity dominates the phase velocity deeply into the SOL. The energy budget is such that RX turbulence pumps expansion free energy from the plasma pressure gradient into flow energy via mode coupling [3,37]- [39]. A similar visualization of density fluctuation structures is given by BES on DIII-D [6] (as shown in figure 6(b)) and GPI on NSTX [7]. From our DIII-D simulations, the poloidal correlation length is ∆ p ∼ 16ρ s ∼ 4 cm, the radial correlation length is ∆ r ∼ 5ρ s ∼ In the far SOL, neutral particle sources make the measured particle flux larger than simulations due to the absence of source particles there in the simulations (reprinted by permission from [41]).
1 cm, and the poloidal propagation velocity is v ph ∼ −10 km s −1 in the edge and v ph ∼ 2 km s −1 in the SOL. The radial correlation length is smaller than the poloidal correlation length by a factor of 4, indicating elongated structures in the poloidal-radial plane as in [6] but with a ratio of 2.5. Figure 7 shows images of the boundary turbulence in the parallel versus perpendicular plane in a magnetic flux surface from BOUT and GPI experiments on NSTX [7]. In the direction parallel to the magnetic field, the correlation length ∆ is very long, ∆ πqR ∼ 10 m. The boundary plasma turbulence thus has a filamentary structure (a set of stripes aligned along field lines) extending only a few cm perpendicular to the magnetic field but many metres in the parallel direction. These filaments as shown in figure 7 from the simulations are qualitatively consistent with the GPI experimental observations. The typical fluctuation parameters obtained are a correlation time τ c ∼ 15-30 µs at the midplane, which is in agreement with probe, BES and the PCI measurements [40], as well as GPI, as shown in figure 8(a). The frequency spectra of the midplane fluctuating density from GPI measurements and BOUT simulations in L-mode in figure 8(b) shows reasonable agreement.
Particle transport perpendicular to the magnetic field Γ r results from correlated fluctuations of the plasma drift velocityṽ E and densityñ, and can be calculated from Γ r = ṽ Eñ . A strong poloidal asymmetry of the turbulent flux of particles is shown in figure 9(b). Due to the X-point null, the plasma drift velocityṽ E ∝ k p ∝ nB t /RB p diverges there as B p → 0, and thus a large particle flux is produced near the X-point region. The particle flux from C-Mod simulation is consistent with Langmuir probe measurement at the midplane near the separatrix as shown in figure 9, with caveat that the probe-inferred flux is a very rough estimate at best due to the perturbation by the presence of the probe body itself.

Summary and conclusions
We have developed a simulation model BOUT for boundary plasma turbulence in toroidal magnetic fusion devices. BOUT contains much of the relevant physics for the boundary problem in the experimentally relevant X-point divertor geometry. Encouraging results have been obtained when comparing with experimental measurements in C-Mod, DIII-D and NSTX. The physics of resistive X-point turbulence and its relation to flow shear generation is discussed. The poloidal fluctuation phase velocity shows structure across the separatrix that is experimentally observed in many fusion devices. The turbulent 3D structures and spectra also resemble the experimental measurements. The particle flux measured from BOUT simulation is consistent with Langmuir probe measurement. A preliminary attempt has been made to characterize motion of 3D turbulent structures.
The BOUT code can have an important impact on the understanding of the dynamics of boundary turbulence. There is strong coupling among the plasma, neutral and impurity species in the SOL. Further development of BOUT including ionization and radiation, combined with more experiments of this nature, will undoubtedly help build the most detailed picture to date of the dynamics of the boundary plasma.